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ABSTRACT 


Low  emissions  of  CO,  NOx,  and  unbumed  hydrocarbons  (UHC)  are  a  difficult  challenge  in  the 
design  of  new  military  gas  turbine  combustors.  Simulation  tools  that  can  predict  emissions  are 
needed  to  reduce  the  cost  of  producing  improved,  low  emissions  combustor  designs.  In  this 
SBIR,  CFD  Research  Corporation  (CFDRC)  continued  to  develop  combustion  Large  Eddy 
Simulation  (LES)  techniques  to  create  a  high  fidelity  tool  for  predicting  emissions.  The  LES 
code  was  improved  by  the  development  and  implementation  of  a  new  multi-step  assumed  PDF 
method  that  accounts  for  more  detailed  kinetics  with  turbulent  chemistry  interactions.  This  new 
method  enables  efficient  turbulent  combustion  CFD  calculations  for  both  steady  state  Reynolds 
Averaged  Navier  Stokes  (RANS)  and  LES  with  multi-step  global  mechanisms.  Tabulation 
methods  were  implemented  and  tested  for  improved  computational  efficiency.  Improvements  to 
the  existing  combustion  models  and  inlet  boundary  conditions  for  LES  were  also  performed.  In 
addition  to  the  new  turbulent  combustion  models,  the  capability  to  generate  the  necessary  global 
mechanisms  from  detailed  reaction  mechanisms  was  developed.  The  final  code  was  validated 
against  benchmark  experimental  data,  and  applied  to  the  Rolls-Royce  JSF  combustor.  Validation 
cases  included  both  premixed  and  diffusion  flames  covering  a  broad  range  of  flame  conditions. 

Although  much  progress  was  made  in  this  Phase  II  effort,  continued  work  is  needed  to  make  the 
new  multi-step  assumed  PDF  model  robust  and  practical.  In  particular,  a  new  solver  for  the 
species  transport  equations  needs  to  be  implemented  to  reduce  run  times  by  a  factor  of  two  or 
more. 


1 


8503/8 


EXECUTIVE  SUMMARY 


Low  emissions  of  pollutants,  including  CO,  NOx,  UHC,  and  smoke,  are  becoming  a  requirement 
for  today's  and  future  military  gas  turbine  engines.  Advanced,  high  performance  gas  turbines 
will  feature  higher  inlet  temperatures  and  pressures,  and  will  require  new  fuel  injector  and 
combustor  designs  if  lower  emissions  are  to  be  realized.  These  low  emissions  gas  turbine 
combustors  will  operate  at  near  stoichiometric  overall  fuel-air  ratios  at  full  power,  and  will  have 
many  emission  design  challenges. 

To  meet  these  challenges  in  a  cost-  and  time-effective  manner,  improved  combustor  design 
systems  are  needed  that  include  accurate  models  of  chemical  kinetics,  turbulent  mixing,  and  their 
interaction.  In  this  SBIR  program,  CFDRC  continued  to  develop  3D  combustion  Large  Eddy 
Simulation  (LES)  capabilities  to  help  meet  these  design  challenges.  Newly  developed 
capabilities  for  capturing  turbulence-chemistry  interactions  with  complex  chemistry,  improved 
multi-step  reactions  mechanisms  for  JP-8,  and  algorithms  to  improve  computational  speed 
greatly  enhance  the  ability  to  predict  emissions  in  high  performance  military  combustors. 

The  primary  focus  of  the  program  was  on  the  development  of  a  new  multi-step,  assumed  PDF 
model  for  turbulent  combustion.  The  newly  developed  model  is  generally  applicable  to  multi- 
step  global  mechanisms  and  results  in  a  minimal  increase  in  the  overall  computational  expense. 
The  effect  of  turbulence  is  accounted  for  by  integrating  the  instantaneous  reaction  rate  over  a 
PDF  to  get  the  resulting  mean  rates.  The  PDF  used  is  based  on  the  mixture  fraction  and  reaction 
progress  variables  to  allow  for  the  characterization  of  both  premixed  and  diffusion  flames.  With 
the  two-variable  PDF  formulation,  oaly  two  additional  transport  equations  must  be  solved  for, 
regardless  of  the  number  of  species  of  steps  in  the  mechanism.  This  capability  represents  a 
significant  improvement  in  the  state-of-the-art  for  practical  turbulent  combustion  modeling. 

Although  the  PDF  model  was  successfully  developed  and  tested,  improvements  in  the  solution  of 
the  transport  equations  will  be  essential  for  the  model  to  find  wide-spread  application.  The 
existing  method  for  solving  the  species  transport  equations  utilizes  a  cell-by-cell  point  solver. 
Specifically,  the  solver  ‘sweeps’  over  all  of  the  computational  cells  in  the  domain  solving  the 
coupled  species  equations  at  each  cell.  While  adequate  for  small  problems,  large  problems 
require  a  large  number  of  iterations  in  order  to  assure  global  convergence.  Steady-state  problems 
can  usually  afford  a  large  number  of  global  iterations  to  achieve  ample  convergence.  However, 
transient  cases  such  as  LES  need  to  keep  the  number  of  global  iteration  per  time  step  to  a 
minimum  if  practical  run  times  are  to  be  realized.  Practical  simulations  will  require  the 
implementation  of  a  “whole-field”  solver  for  the  coupled  species  transport  equations. 

In  addition  to  the  new  turbulent  combustion  models,  the  ability  to  generated  global  mechanisms 
from  detailed  reaction  mechanisms  was  also  developed.  Multi-step  kinetics  are  essential  to 
accurately  capture  behavior  of  real  reactions,  such  as  the  conversion  of  CH4  to  CO  and  then  to 
CO2.  The  ability  to  generate  these  mechanisms,  commonly  termed  Chemical  Reactor  Modeling 
(CRM),  was  shown  to  an  effective  method  of  determining  global  Arrhenius  rate  parameters.  The 
global  mechanisms  generated  can  be  fine  tuned  to  combustor  operating  conditions  including  inlet 
temperature,  equivalence  ratio  range,  and  operating  pressure. 
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1. 


INTRODUCTION 


Low  emissions  of  pollutants,  including  CO,  NOx,  UHC,  and  smoke,  are  becoming  a  requirement 
for  today’s  and  future  military  gas  turbine  engines.  Advanced,  high  performance  gas  turbines 
will  feature  higher  inlet  temperatures  and  pressures,  and  will  require  new  fuel  injector  and 
combustor  designs  if  lower  emissions  are  to  be  realized.  These  low  emissions  gas  turbine 
combustors  will  operate  at  near  stoichiometric  overall  fuel-air  ratios  at  full  power,  and  will  have 
many  emission  design  challenges,  some  of  which  are: 

1 .  How  to  effectively  dilute  fuel-rich,  high  CO  regions  before  the  hot  gas  enters  the  turbine; 

2.  How  to  minimize  residence  times  for  high-temperature,  thermal  NOx  producing  regions; 
and, 

3.  How  to  oxidize  CO  and  reduce  UHC  at  low  power  conditions. 

To  meet  these  challenges  in  a  cost-  and  time-effective  manner,  improved  combustor  design 
systems  are  needed  that  include  accurate  models  of  chemical  kinetics,  turbulent  mixing,  and  their 
interaction.  Current  steady-state  RANS  (Reynolds  Averaged  Navier  Stokes)  CFD  combustion 
codes  lack  the  quantitative  accuracy  needed  for  reliable  predictions  of  emissions.  The  limitations 
of  the  current  CFD  codes  are  mainly  due  to  deficiencies  in  the  treatment  of  the  turbulent  fluid 
motion  and  its  interaction  with  chemical  kinetics,  i.e.,  turbulence-combustion  interaction.  This  is 
because  turbulent  combustion  is  inherently  unsteady  (including  vortex  shedding,  shear  layer 
mixing,  and  acoustic  wave  propagation),  and  RANS  codes  cannot  capture  counter-gradient 
diffusion  and  other  unsteady  phenomena  seen  in  gas  turbine  combustors.  The  accuracy  of 
emissions  predictions  is  especially  vulnerable  to  RANS  limitations.  In  this  SBER  program, 
CFDRC  is  continuing  to  develop  3D  combustion  LES  capabilities  to  help  meet  these  design 
challenges.  Newly  developed  capabilities  for  capturing  turbulence-chemistry  interactions  with 
complex  chemistry,  improved  multi-step  reactions  mechanisms  for  JP-8,  and  algorithms  to 
improve  computational  speed  will  greatly  enhance  the  ability  to  predict  emissions  in  high 
performance  military  combustors,  including  the  accurate  prediction  of  NOx,  CO,  and  UHC. 

In  Phase  I,  the  feasibility  of  improved  emissions  predictions  for  CO  and  NOx  using  combustion 
LES  modeling  of  an  aero  gas  turbine  combustor  was  demonstrated.  The  combustion  LES 
analysis  utilized  a  local  subgrid  turbulent  kinetic  energy  model  (called  the  localized  dynamic 
kinetic  energy  model,  LDKM)  to  predict  local  viscosities  as  a  function  of  the  grid  size  and 
velocity  gradients.  The  turbulent  kinetic  energy  provided  a  source  term  for  the  subgrid 
turbulence-chemistry  assumed  PDF  model.  LES  filtered  reaction  rates  for  the  fuel,  CO,  and  NOx 
transport  equations  included  turbulent  fluctuations  in  the  mixture  fraction  and  reaction  progress 
variables.  A  2-step  chemical  reaction  model  that  provides  a  finite  rate  of  fuel  oxidation  to  CO 
and  then  CO  oxidation  to  equilibrium  products  was  used.  Also,  a  de-coupled  NOx  transport 
method,  accounting  for  thermal  NOx  reactions,  was  included. 

Both  combustion  LES  and  steady-state  RANS  calculations  of  the  Rolls-Royce  AE3007 
combustor  were  performed.  The  2.3  million  cell  computational  domain  consisted  of  the  entire 
combustor  internal  flowfield  (from  the  swirl  vane  exit  planes  to  the  combustor  exit),  as  well  as 
the  outer  annuli  surrounding  the  combustor.  Estimated  flowsplits  were  provided  by  Rolls-Royce. 
Comparisons  of  predictions  to  available  experimental  data  at  the  combustor  exit  plane  were 
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made.  Figures  1  and  2  show  snapshots  and  time-averaged  predictions  of  velocity  magnitude  and 
temperature  (high-power  conditions).  The  velocity  predictions  indicate  stronger  penetration  of 
the  primary  jets  using  LES.  This  stronger  penetration  closes  off  the  central  recirculation  zone, 
making  it  more  compact  for  the  LES  predictions.  The  temperature  predictions  show  the  effect  of 
stronger  mixing  with  LES.  The  peak  temperature  is  reduced  using  LES  and  the  hot  temperatures 
spread  out  over  a  greater  area  in  the  intermediate  zone.  This  higher  mixing  in  the  intermediate 
zone  is  due  to  the  strong  primary  jets  that  create  stronger  wakes  and  mixing  behind  the  jets. 


LES  Snapshot 


LES  Time-Average  LES  Time-Average 


Figure  2.  LES  and  RANS  Predictions  of 
Temperature  and  Spray  Trajectories  at  High 
Power  Conditions 


LES  Snapshot 


(m/») 


Figure  1.  LES  and  RANS  Predictions  of 
Velocity  Magnitude  at  High  Power  Conditions 
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Figure  3  shows  the  NOx  emissions  predictions  at  high  power  where  most  NOx  is  formed.  NOx  is 
formed  at  high  temperature  locations  when  there  is  excessive  oxygen,  typically  at  local 
equivalence  ratios  around  0.8  -  0.9.  The  NOx  emissions  results  show  that  NOx  is  formed  sooner 
and  over  a  larger  volume  for  the  LES,  particularly  out  towards  the  inner  and  outer  liner  of  the 
intermediate  zone.  However,  higher  peak  values  are  predicted  for  the  RANS  along  the  centerline 
of  the  dilution  zone,  most  likely  due  to  the  higher  peak  temperatures  in  this  region.  Overall, 
more  NOx  is  predicted  for  RANS  at  the  combustor  exit. 

Table  1  shows  comparisons  of  the  predictions  and  measurements  at  the  high  power  conditions. 
Overall,  LES  provides  a  better  prediction  of  emissions  compared  to  the  RANS.  The  improved 
results  can  be  largely  attributed  to  the  stronger  mixing  of  primary  and  dilution  air  with  the 
combustor  core  flow  for  the  LES  and  the  mixing  produced  by  combustor  unsteadiness.  The 
RANS  did  not  capture  the  higher  jet  penetration  as  observed  in  the  LES.  The  CO  emissions 
show  the  largest  differences  between  the  two  predictions  and  measurements.  However,  the 
absolute  magnitude  of  the  CO  is  small,  making  quantitative  predictions  difficult. 


LES  Snapshot  LES  Time- Average 


Figure  3.  Comparison  ofNOx  Predictions  for  LES  and  RANS  at  High  Power  Conditions 
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Table  1.  Predicted  and  Measured  CO,  NOx,  and  Temperature  for  AE3007  High  Power 

Conditions 


Experimental 

RANS 

LES 

T/Tmeasured 

1 

0.98 

1.00 

CO/COmeasured 

1 

10.73 

0.42 

N Ox/N  0  x f  measured 

1 

1.31 

1.15 

Computations  were  next  carried  out  for  low  power  conditions.  The  CO  emissions  predictions  for 
this  calculation  are  shown  in  Figure  4.  The  CO  emissions  at  low  power  conditions  can  be 
significant,  since  quenching  for  CO  oxidation  reactions  can  occur  at  the  lower  combustor 
temperatures.  Rolls-Royce  combustor  exit  measurements  indicated  ~30  times  more  CO  at  low 
power  compared  to  high  power  conditions,  but  still  at  relatively  low  levels.  In  the  CFD 
predictions,  a  large  amount  of  CO  is  predicted  in  the  primary  zone,  and  some  oxidation  occurs  in 
the  intermediate  zone.  However,  complete  oxidation  does  not  occur  in  the  intermediate  zone, 
and  the  dilution  air  quenches  the  CO  at  nonequilibrium  values  in  the  dilution  zone.  The  expected 
cause  of  high  CO  at  lower  power,  i.e.,  liner  cooling  quenching,  is  not  evident  in  these 
calculations,  probably  because  the  high  CO  in  the  mainstream  masks  any  quenching  that  may 
occur  along  the  liner.  The  CO  predictions  results  seem  to  indicate  that  the  LES  with  stronger 
mixing  reduces  the  CO  emissions  at  the  combustor  exit  compared  to  RANS.  This  is  most  likely 
due  to  the  shortened  slightly  fuel-rich  primary  zone  predicted  in  the  LES  case.  The  exit  CO 
emissions  for  LES  are  -20%  less  than  the  RANS  predictions.  Table  2  shows  comparisons  of  the 
predictions  with  experimental  data.  Despite  the  improvement  for  LES,  both  RANS  and  LES 
over  predict  the  combustor  exit  CO  emissions  by  substantial  amounts.  It  is  expected  that 
improvements  to  the  chemical  kinetics  mechanism  will  help  the  low  power  CO  predictions.  The 
current  CO  oxidation  model  assumes  equilibrium  levels  of  OH  and  does  not  include  finite-rate 
kinetics  of  other  minor  species. 
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LES  Snapshot 


LES  Time-Average 


Figure  4.  Comparison  of  CO  Emissions  Predictions  for  LES  and  RANS  at  Low  Power 

Conditions 


Table  2.  Predicted  and  Measured  CO,  NOx,  and  Temperature  for  AE3007  Low  Power 

Conditions 


Experimental 

RANS 

LES 

T/TmeasUred 

1 

0.94 

0.96 

C  O/C  O  mea5  ured 

1 

9.27 

7.5 

NOjt/NOx,  measured 

1 

1.12 

0.85 

The  difference  between  the  LES  and  RANS  NOx  emissions  at  low  power  was  relatively  small.  It 
can  be  seen  that  the  NOx  is  less  at  the  combustor  exit  for  the  LES.  The  lower  predicted  NOx  for 
LES  is  due  to  the  better  mixing  and  reduction  in  peak  temperature  for  the  LES.  The  comparison 
with  data  in  Table  2  indicates  that  the  RANS  is  -12%  higher  than  data  and  the  LES  is  -15% 
lower  than  the  measurements.  Since  the  measured  NOx  is  quite  low,  these  small  discrepancies 
are  not  too  serious  and  both  LES  and  RANS  essentially  provide  equally  good  predictive 
capabilities  for  low  power  NOx. 
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The  Phase  I  research  demonstrated  the  combustion  LES  has  significant  potential  to  provide 
improved  capabilities  to  accurately  predict  emissions.  Needed  capabilities  with  respect  to 
combustion  modeling  were  identified.  Specific  areas  identified  during  Phase  I  where 
improvements  in  modeling  capability  are  needed  include  turbulence  chemistry  interactions  with 
multi-step  chemical  mechanisms,  relevant  multi-step  mechanisms  for  JP-8,  and  algorithms  to 
improve  the  speed  of  the  complex  chemistry  calculations.  Work  performed  in  Phase  II  addresses 
needs  related  to  turbulence-chemistry  interaction  and  multi-step  chemistry.  Tasks  include 
development  and  validation  of  the  new  capabilities,  along  with  application  to  the  Rolls  Royce 
JSF  combustor. 

2.  TECHNICAL  OBJECTIVES 

The  overall  goal  of  this  SBIR  project  was  to  develop  an  analysis  tool  that  captures  the 
fundamental  physics  of  NOx,  CO,  and  UHC  emissions  in  high  performance  gas  turbine 
combustors.  In  Phase  I,  the  feasibility  of  the  proposed  analysis  to  predict  emissions  was 
demonstrated  by  predicting  CO  and  NO*  emissions  in  a  production  Rolls-Royce  aeroengine  gas 
turbine  combustor.  In  Phase  II,  the  development  of  the  LES  combustion  models  continued.  A 
significant  amount  of  validation  took  place  and  the  capability  was  applied  to  the  JSF  combustor. 
Specific  objectives  of  the  current  Phase  II  program  were: 

1.  Incorporate  multi-step  assumed  PDF  method  into  combustion  LES  code. 

2.  Develop  and  apply  reduced  (10-12  species)  JP-8  chemical  kinetic  mechanism  for  JSF 
combustor  conditions. 

3.  Validate  the  combustion  LES  code  against  benchmark  DOE  SimVal  combustor  data  for 
predicting  combustion  dynamics  and  emissions  at  lean  conditions  using  gaseous  fuels. 

4.  Validate  the  combustion  LES  code  against  published  validation  data  for  premixed  and 
diffusion  flames. 

5.  Apply  combustion  LES  code  to  the  Rolls-Royce  JSF  combustor  design. 

Objectives  one  through  five  were  achieved  during  the  completion  of  the  program.  However,  the 
validation  and  application  of  the  combustion  model  were  more  limited  than  originally  planned, 
due  to  difficulties  in  solving  the  coupled  enthalpy  and  species  equations.  The  technical 
challenge  associated  with  the  solution  procedure  was  discovered  late  in  the  Phase  II  effort  and 
the  solution  of  the  problem  was  beyond  the  scope  of  the  Phase  II  effort.  Sufficient  simulations 
were  performed  to  demonstrate  the  feasibility  of  the  multi-step  PDF  combustion  model,  but 
application  of  the  capability  to  transient  LES  simulations  will  be  limited  until  the  inadequacies  of 
the  solution  procedure  is  addressed. 

3.  RESULTS 

The  following  section  contains  a  complete  description  of  the  work  performed  during  this  Phase  II 
SBIR  project.  Table  3  shows  a  graph  of  the  project  schedule  and  marks  the  completion  of  the 
project  milestones. 
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Table  3.  Work  Schedule 


2  Perform  3D  Calculation  of  Rolls-Royce  JSF  7  Complete  Second  Validation  Case  I  I  Planned 

Combustor  Using  Existing  Capabilities  8  Complete  Model  Improvements  _ 

3  Complete  ISAT  Implementation  with  Multi-Step  PDF  9  Complete  Baseline  Simulation  of  JSF  Combustors  - — *  Perform* 

4  Complete  JP8  Reduced  Mechanism  10  Complete  Analysis  of  JSF  Combustors 

5  Complete  FirstValidation  Case 


Extensive  improvements  in  the  combustion  modeling  capabilities  were  performed  in  CFD-ACE+ 
during  this  Phase  II  SBIR  program.  All  the  efforts  were  focused  on  capabilities  applicable  to 
Large  Eddy  Simulation  (LES),  or  that  can  be  used  for  both  RANS  and  LES  calculations.  Section 
3.1  describes  the  improvements  and  models  added  to  the  existing  LES  capabilities,  and  also 
provides  an  overview  the  LES  modeling  capabilities  in  CFD-ACE+.  The  new  multi-step  PDF 
combustion  model  developed  under  this  SBIR  program  is  described  in  Sections  3.2. 
Improvements  in  LES  boundary  conditions  are  discussed  in  Section  3.3.  Generation  of  multi- 
step  mechanisms  is  discussed  in  Section  3.4,  followed  by  the  discussion  of  the  validation  results 
in  Section  3.5.  The  results  of  the  JSF  simulations  are  briefly  described  in  Section  3.6,  and  the 
results  of  these  analyses  are  included  in  an  export-controlled  and  proprietary  addendum  report. 


3.1  Improvements  to  Existing  LES  Models 

LES  Overview:  LES  is  a  time-accurate  solution  to  the  Navier  Stokes  equations  with  grid  sizes 
and  time  steps  that  are  small  enough  to  resolve  the  energy-containing  scales  of  a  turbulent  flow. 
Subgrid  models  for  turbulent  mixing  and  chemical  reaction  are  needed  to  resolve  the  small-scale 
effects  that  are  more  universal  and  easier  to  approximate.  Each  variable  in  a  LES  is  decomposed 
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into  a  large-scale  component  (indicated  by  an  overbar)  and  a  residual  component  (indicated  by  a 
prime), 


f  =  f +  f' 


The  large-scale  component  is  obtained  by  spatially  averaging  f  with  a  filter  function  G, 

f  =  Jf(x')G(x,x')dx' 


(1) 

(2) 


The  filter  effectively  eliminates  fluctuations  on  scales  smaller  than  a  specified  size.  The  large 
scale,  or  resolved  components  are  time  dependent,  in  contrast  to  the  average  components  in 
conventional  (Reynolds  averaged)  turbulence  modeling.  For  variable  density  flows,  CFDRC 
uses  Favre  (density  weighted)  filtering: 


f 


(3) 


With  the  application  of  an  LES  Favre  averaging  procedure  to  the  governing  transport  equations 
(mass,  momentum,  energy,  species  mass  fractions,  respectively),  we  obtain: 


dp  ,  dpu;  _Q 
<3t  dxj 

5puj  SpUjU  j  oTjj  STjj 

dt  dxj  dxj 

dpE  dpujE  _  dijjUj 

dt  3xj  3xj  5x, 


(4) 


dp^q  dpUj$a  dJ”  dM“ 

dt  5x[  9xj  5xj  a 

The  filtered  equations  contain  unknown  terms  such  as  t;j  (velocity-velocity  correlation)  arising 
from  the  filtering  of  nonlinear  terms  and  are  known  as  subgrid  scale  (SGS)  stresses.  The  closure 
problem  of  turbulent  flows  using  LES  is  handled  by  providing  models  for  the  various  subgrid 
scale  correlations. 

Subgrid  Turbulence:  Unlike  typical  steady-state  turbulence  models,  the  turbulence  models  for 
LES  compute  an  eddy  viscosity  that  is  a  function  of  the  grid  (or  filter)  size.  The  larger  the  size 
of  the  computational  cells  the  higher  the  value  of  the  subgrid  turbulent  viscosity.  Thus,  as  the 
grid  is  made  finer,  the  modeled  effect  of  subgrid  turbulent  mixing  becomes  less  important  and 
more  large-scale  phenomena  are  then  directly  computed.  The  Smagorinsky  model 
(Smagorinsky,  1963)  computes  the  turbulent  viscosity  from  the  magnitude  of  the  resolved  strain 
tensor  Sy,  the  grid  filter  width  A,  and  the  Smagorinsky  constant  Cs  according  to: 
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(5) 


vt  =(CSA)2|S|,  whereSij 


f  duj 

dx : 
J 


+ 


du£ 
ax;  J 


where  the  grid  filter  width  is  computed  as  the  cube  root  of  the  grid  cell  volume.  This 
implementation  allows  LES  calculations  with  arbitrary  grid  types,  including  hybrid  grid 
schemes.  Various  studies  have  shown  that  the  tunable  parameter  Cs  in  the  Smagorinsky  model  is 
not  a  universal  constant  (Avva  and  Sundaram,  1998).  It  has  been  suggested  that  the  Cs  parameter 
should  be  flow  dependent  and  should  vary  from  region  to  region  in  complex  flows. 

A  LES  model  which  provides  a  flow  dependent  subgrid  viscosity  parameter  is  Localized 
Dynamic  subgrid  Kinetic  energy  model  (LDKM)  and  is  the  basis  for  the  newly  developed  model 
for  the  scalar  dissipation  coefficient.  The  LDKM  model  was  developed  by  Kim  and  Menon 
(1997)  to  provide  subgrid  stresses  without  a  priori  specification  of  any  constants.  The  LDKM 
uses  scale  similarity  and  the  subgrid-scale  kinetic  energy 


k 


sgs  - 


2 


(ukuk  -ukuk) 


(6) 


to  model  the  unresolved  scales.  Using  ksgs  the  subgrid  stress  tensor  is  modeled  as 


r - - OC  Ak'^S-  +-8  -lc 

MJ  ~  z'v-'T£mSgS  °1J  ^  ^  °lJASgS 


(7) 


with  the  resolved-scale  strain  tensor  defined  as 


sij  =- 
1J  2 


(Tu;  <9u; 
'-  +  ■ — - 


dx  j  3xj 


(8) 


In  the  modeling  of  the  SGS  stresses,  implicitly  the  eddy  viscosity  is  parameterized  as 


vT  =CTAkJ2  • 


(9) 


The  subgrid-scale  kinetic  energy  is  obtained  by  solving  the  transport  equation 


^sgs  _  ^sgs  3u:  d 

dt  dx-,  J  dx :  8  dx-. 


lJ 


vt 


i  v 


dk sgs^ 

5Xj  j 


GO) 


which  is  closed  by  providing  a  model  for  the  subgrid  dissipation  rate  term,  ssgs,  based  on  simple 
scaling  arguments 
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(11) 


In  these  models,  CT  and  CE  are  adjustable  coefficients  determined  dynamically  using  the 
information  from  a  resolved  test-scale  field.  The  test-scale  field  is  constructed  from  the  large- 
scale  field  by  applying  a  test  filter,  which  is  characterized  by  A ,  the  test  filter  width.  In  the  LES 
code,  with  arbitrary  grids,  we  are  using  a  test  filter  consisting  of  a  weighted  average  of  the  cells 
sharing  a  node  with  the  current  cell.  This  average  is  biased  towards  the  current  cell,  with  a 
weight  equal  to  the  number  of  vertices  of  the  cell.  The  cells  that  share  a  face  with  a  current  cell 
have  a  weight  of  two. 

The  application  of  the  test  filter  on  any  variable  is  denoted  by  the  top  hat.  By  definition,  the 
Leonard  stress  tensor  at  the  test-scale  level  is 

Lij  =u[uj-UiUj  (12) 

The  Leonard  stress  tensor  and  the  SGS  tensor  are  known  to  have  high  degrees  of  correlation, 
which  justifies  the  use  of  similarity  in  the  derivation  of  the  dynamic  model  coefficients.  The 
resolved  kinetic  energy  at  the  test  filter  level  is  defined  from  the  trace  of  the  Leonard  stress 
tensor 


ktest  =^(ukUk-ukuk)  ■ 


This  test  scale  kinetic  energy  is  dissipated  at  small  scales  by 


(13) 


e  =  (v  +  vx) 


eft,  <3u j  c'uj  c'uj 
<3xj  dxj  dxj  dxj 


(14) 


Based  on  a  similarity  assumption  and  using  appropriately  defined  parameters,  the  Leonard  stress 
tensor  has  a  representation  analogous  to  the  SGS  stress  tensor 


Lij  =-2CtAk^,Sjj  +lf>ijL|tjc 


(15) 


The  least  square  method  is  applied  to  obtain  the  model  constant 

1  Lijaij 


CT  = 


2  CTjjCTjj 


(16) 


where 
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(17) 


o-=-Akl/2S- 

aU  ^test^ij 


Finally,  a  corresponding  approach  is  used  to  determine  the  dissipation  rate  constant.  By  invoking 
similarity  between  the  dissipation  rates  at  the  subgrid  level  and  at  the  test  scale  level  C6  is 
determined  to  be 


(dUj/dXj)(dUi/dXj)- 

[(v  +  vT  )a]~'  k^t 


(18) 


The  coefficients  of  the  LDKM  model  are  Galilean  invariable  and  realizable.  This  model  is  also 
quite  simple  and  efficient,  does  not  rely  on  any  ad  hoc  procedures,  and  it  is  applicable  to  various 
flow  fields  without  adjustment  of  the  model.  The  subgrid  turbulence  kinetic  energy  ksgs  from 

the  LDKM  is  the  required  input  for  several  subgrid  chemistry  models  such  as  the  assumed  PDF 
and  the  Linear  Eddy  Model. 

3.1.2  Development  of  LDVM  for  the  PDF  Combustion  Model 

The  equations  solved  by  CFD-ACE+  for  LES  of  turbulent  reacting  flows  are  transport  equations 
for  filtered  density-weighted  values.  The  effect  of  turbulence  on  chemical  reaction  and  on 
composition  dependent  variables,  such  as  density  or  temperature,  must  be  accounted  for.  The 
filtered  values  of  density  and  temperature  cannot  be  calculated  from  the  filtered  value  of  the 
mass  fractions  since  they  are  nonlinear  functions  of  the  composition.  These  variables  can  be 
calculated  more  accurately  by  integrating  the  product  of  the  variable  of  interest  and  the  filtered 
density-weighted  joint  composition  probability  density  function  (PDF)  over  the  range  of 
composition  values.  In  order  to  have  a  simpler  model  to  work  with  and  to  increase  computational 
efficiency,  we  have  decided  to  focus  on  the  assumed  PDF  approach. 

The  filtered  joint  composition  PDF  is  a  complete  statistical  description  of  the  composition  of  the 
fluid  at  a  single  point  in  space  and  time.  If  the  PDF  is  known,  then  the  filtered  value  of  any 
function  of  composition  can  be  evaluated  by  multiplying  that  function  by  the  PDF  and 
integrating  over  the  range  of  possible  compositions. 

co(Y1,...,YN;x,t)  =  fco(Y1,...,YN) 

(P(Y1,...,YN;x,t)dYl...dYN)  (19) 


where  P(Yi,..,YN|x,t)  is  the  filtered  joint  PDF  of  the  N  mass  fractions  at  the  position  x  and  time  t 
and  co(Yi,...,Yn)  is  an  arbitrary  function  of  the  mass  fractions.  Favre-filtered  quantities  can  be 
calculated  by  defining  a  filtered  density-weighted  PDF. 


P(Y1,...,Yn) 


Ip(Yi„..,Yn)p(Yi . yn) 

p(Y,,.  ,Yn) 


(20) 
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A  parametric  form  of  the  PDF  is  assumed  and  the  parameters  in  the  model  are  related  to 
variables  governed  by  transport  equations.  The  parametric  form  of  the  PDF  assumes  the 
composition  can  be  specified  by  a  single  mixture  fraction  and  a  single  reaction  progress  variable. 
This  assumption  limits  the  reaction  models  available  when  the  prescribed  PDF  model  is  used. 
The  mass  diffusivities  of  all  species  must  be  equal  and  no  more  than  two  mixtures  can  be 
defined.  In  the  following  sections,  determination  of  the  PDF  and  the  filtered  variables  will  be 
described. 

Determination  of  PDF:  The  joint  composition  PDF  is  a  function  of  a  mixture  fraction  and  a 
reaction  progress  variable.  The  reaction  progress  is  defined  as 

9  (Yf.mtt.-Yf) 

(Yf.ttax-Yf.min) 

where  Yf  is  the  mass  fraction  of  the  fuel  in  the  one  step  reaction  and  the  minimum  and 
maximum  values  are  functions  of  the  mixture  fraction.  The  mixture  fraction  and  reaction 
progress  are  assumed  to  be  independent,  so  the  two-dimensional  PDF  is  a  product  of  the  two 
one-dimensional  PDF's. 


P(f,0)  =  P(f)P(0) 


(22) 


The  one-dimensional  PDF's  both  have  two  parameters  that  are  related  to  the  filtered  mean  and 
the  subgrid  variance  of  the  mixture  fraction  or  reaction  progress. 

In  the  present  approach,  transport  equations  are  solved  for  both  the  filtered  mean  and  variance  of 
the  corresponding  variable.  The  transport  equations  for  the  variances  of  the  mixture  fraction  and 
reaction  progress  include  production  terms  caused  by  gradients  in  the  filtered  values,  dissipation 
terms,  and  (for  the  reaction  progress)  a  term  due  to  chemical  reaction. 


3-  2 

apCf  + 


—  2  3 

pUjaf=- 


ao?l 

SxJ 


+  2T 


df 


dx 


r  p8~2 


J  J 
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—  r  P£  o 

4  k  a0+ 


(23) 


(24) 


In  the  preceding  equations,  Cx  has  the  value  of  2  in  RANS  methodology.  To  improve  the 
predictions  in  LES  we  have  derived  a  Localized  Dynamic  Subgrid  Variance  model  (LDVM). 
This  new  closure  follows  the  logic  behind  the  derivation  the  Localized  Dynamic  Subgrid  Kinetic 
energy  model  (LDKM).  Furthermore,  computationally  the  new  subgrid  variance  model  relies  on 
information  provided  by  the  LDKM  with  quantities  such  as  subgrid  kinetic  energy  and  subgrid 
kinetic  energy  dissipation. 
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In  the  new  model,  Cx  is  an  adjustable  coefficient  determined  dynamically  using  the  information 
from  a  resolved  test-scale  field.  The  derivation  starts  from  the  definition  of  the  variance  of  scalar 
Z  (representing  either  the  mixture  fraction  or  the  reaction  progress  variable).  The  subgrid  scalar 

variance  is  defined  as  ViSs  =  Z2  -Z  ,  whereas  the  scalar  variance  at  the  test  level  is 

Vtest  =  Z  -  Z  .  From  dimensional  analysis  the  subgrid  scalar  dissipation  is  modeled  by 


Xsgs  s  2D 


'dZ_  dZ 
KdXj  dxt 


dZ  dZ  ' 
dxi  dxi  j 


(25) 


and  by  invoking  the  principle  of  similarity  between  the  test  level  and  the  subgrid  level  quantities, 
the  scalar  dissipation  at  the  test  level  is 


f 


Xtest  —  +  Dj  ) 

v 


&Z  dZ 
dxj  5xj 


az  az^ 

5x;  dx: 

1  V 


(26) 


where  DT  is  the  turbulent  diffusivity.  Since  xmt  is  fully  resolved  at  every  point  in  the 
computational  domain,  it  will  be  used  to  calculate  the  unknown  parameter  C  .  The  turbulent 
dissipation  at  the  test  level  is  expressed  as 


8  test 


(27) 


Using  equation  (27),  from  equation  (26)  the  final  expression  for  the  dynamic  coefficient  Cz  is 
obtained  as 


2(D  +  Dy  )A 


CX 


r  &L  cfL 
5x;  dx. 


^  \ 

dZ  dZ 
dxj  dxj 


test^test 


(28) 


This  generic  expression  is  used  in  the  transport  equations  of  the  mixture  fraction  variance  and 
reaction  progress  variance  to  compute  the  model  coefficient  of  the  dissipation  term.  In  the  next 
section  we  detail  the  results  obtained  with  this  new  model  for  a  bluff-body  combustor. 

3.1.3  Testing  of  LDVM  Model 

The  LES  calculations  shown  here  have  utilized  the  localized  dynamic  subgrid  kinetic  energy 
model  (LDKM)  for  subgrid  turbulence  with  1-step  chemistry  and  assumed  PDF  with  LDVM 
model  for  subgrid  chemistry.  Figure  5  shows  a  comparison  between  the  predictions  of  the 
assumed  PDF  model  without  (i.e.  by  using  a  constant  Cz=2)  and  with  LDVM,  for  instantaneous 
temperature  and  the  reaction  progress  variance. 
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Instantaneous  Temperature 
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Progress  Variable  Variance 


Progress  Variable  Variance  (LDVM) 

Figure  5.  Predictions  of  Temperature  and  Variance  in  a  Lean  Premixed  Bluff-Body  Combustor 

Using  Assumed  PDF  Model  with  LES 
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The  most  striking  difference  between  the  two  models  is  apparent  in  the  progress  variable 
variance  field  snapshots.  The  LDVM  model  is  tracking  the  flame  front  very  closely,  which  is  the 
ideal  behavior  for  such  chemistry  models.  Peak  magnitudes  are  also  similar  to  those  seen  in  tire 
flame  location  between  the  two  simulations.  The  2D  calculations  demonstrated  that  the  LDVM 
technique  was  providing  reasonable  results,  both  in  terms  of  location  and  magnitude. 

The  LDVM  model  has  also  been  applied  to  the  3D  DOE  SimVal  combustor.  Work  with  the 
SimVal  combustor  at  DOE  is  designed  to  provide  experimental  data  that  can  be  used  to  validate 
combustion  CFD  codes,  with  particular  emphasis  on  understanding  combustion  instability, 
emissions,  and  variable  fuel  effects  at  actual  gas  turbine  combustor  conditions.  Detailed 
experimental  data  are  being  obtained,  including  emission  images,  velocity,  temperature,  and 
species  maps,  and  dynamic  wall  pressure. 

The  baseline  geometry  of  the  SimVal  combustor  includes  well-defined  acoustic  boundaries  with 
a  choke  plate  immediately  upstream  of  the  swirl  vanes  and  a  choked  nozzle  at  the  downstream 
end  of  a  resonant  section.  Figure  6  shows  the  baseline  geometry.  The  corresponding  CFD 
model  is  shown  in  Figure  7.  The  baseline  operating  conditions  included: 

Air  Mass  Flow-Rate  0.26  kg/sec 
Inlet  Air  Temperature  600  K  (620  F) 

Equivalence  Ratio  0.47  -  0.7 

Pressure  ~  4.5  atm  (varied  with  equivalence  ratio) 
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Figure  6.  Baseline  Geometry  for  SimVal  Combustor 
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Figure  8  shows  the  preliminary  results  of  the  vorticity  magnitude  and  NOx  source  term  at  an 
equivalence  ratio  of  0.7,  obtained  prior  to  the  development  of  the  LDVM  model.  The  NOx  rate 
increases  in  the  shear  layer  where  vorticity  is  high  and  then  decreases  near  the  cool  walls. 


SimVal  Combustor;  PHI=0.7 


SimVal  Combustor;  PHI=0.7 


<T\|  (O  CO  O  CN 


Figure  8.  3D  LES  Snapshot  of  Filtered  NOx  Source  Term  and  Vorticity  Magnitude 

(t=0.7) 


The  NOx  and  CO  emissions  were  over  predicted,  possibly  due  to  improper  estimated  wall 
temperatures.  At  higher  fuel-air  ratios  the  LES  predictions  also  predicted  flashback  into  the 
premix  barrel,  while  the  experiments  did  not. 

Simulations  of  the  SimVal  combustor  were  completed  using  the  LDVM  model  with  the  two- 
variable  PDF  (single-step  rate).  Figure  9  and  10  show  the  progress  and  progress  variance  in  the 
combustor.  The  progress  variable  tracks  the  flame  location  in  the  premixed  combustor.  Both 
plots  show  reasonable  characteristics  using  the  LDVM  model.  Also,  the  tendency  of  the  flame  to 
propagate  along  the  inside  diameter  of  the  premix  barrel  has  been  eliminated.  The  flame  now 
stabilizes  at  a  location  just  inside  in  the  barrel.  This  is  shown  in  Figure  1 1  where  temperature  is 
plotted  on  an  isosurface  of  progress  (0.9  -  corresponding  to  nearly  complete  reaction). 
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Figure  9.  Reaction  Progress  with  the  LDVM  Model 


Figure  10.  Reaction  Progress  Variance  with  the  LDVM  Model 
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Figure  11.  Temperature  Plotted  at  a  Reaction  Progress  of  0.9 


Snapshots  of  temperature  contours  at  various  times  using  the  LDVM  model  are  shown  in  Figure 
12.  Similar  to  the  results  shown  in  Figure  11,  it  can  be  seen  that  the  flame  anchors  on  the 
centerbody,  a  short  distance  upstream  of  the  exit.  This  location  of  flame  anchoring  closely 
mimics  the  experiment  based  on  observations  of  the  experimentalists.  Another  feature  seen  in 
Figure  11  is  the  fine-scale  structure  of  the  flame,  which  is  more  in-line  with  what  is  expected  in 
LES  calculations.  This  is  the  result  of  inlet  turbulence,  which  was  not  captured  in  previous  LES 
runs.  The  improved  boundary  conditions  (described  in  Section  3.3)  effectively  introduce 
reasonable  length  scales  in  the  turbulent  flow  structure  that  propagate  through  the  premix  tube. 

Snapshots  of  NOx  contours  at  various  times  are  shown  in  Figure  13.  NOx  at  the  exit  is 
approximately  80  ppm,  much  higher  than  extrapolated  from  the  measurements.  Unfortunately, 
the  wall  heat  losses  were  not  correctly  accounted  for,  so  the  LES  calculation  had  higher  bulk 
temperatures  in  the  combustor  and  exhaust  duct  than  the  experiment.  Capturing  heat  loss  with 
LES  is  highly  dependant  on  the  wall  function  models  and  near  wall  grid  distribution.  In  an 
ongoing  SBER  program  with  the  Navy,  CFDRC  is  improving  near  wall  heat  transfer  and  velocity 
models  for  LES  simulations  to  provide  accurate  heat  transfer  characteristics. 

The  combustor  computational  model  was  very  quiet  in  terms  of  pressure  dynamics.  Pressure  rms 
values  of  ~0.3  psi  were  predicted.  This  was  in-line  with  extrapolated  pressure  rms 
measurements. 
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Figure  12.  Snapshots  of  Temperature  Contours  in  the  SimVal  Combustor 


NOx  (ppm) 


Figure  13.  Snapshots  ofNOx  Contours  in  the  SimVal  Combustor 
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3.2  Multi-Step  PDF  Turbulent  Combustion  Model 


Particularly  for  gas  turbine  combustion  modeling,  there  is  a  pressing  need  for  a  practical, 
computationally  efficient  turbulent  combustion  model  that  is  applicable  to  all  types  of  flames, 
accounts  for  finite  rate  effects,  and  can  be  used  with  multi-step  chemical  mechanisms.  To  meet 
this  need,  a  new  assumed,  two-variable,  joint-PDF  combustion  model  for  multi-step  reaction 
mechanisms  has  been  developed  in  this  SBIR.  The  model  allows  detailed  chemistry  for  species 
and  combustion  predictions,  with  turbulence-chemistry  interaction  modeled  by  a  joint  PDF  of 
two  independent  variables:  the  mixture  fraction  and  reaction  progress.  A  similar  approach  has 
been  successfully  used  by  Bohn  and  Lepers  (1999).  However,  they  used  temperature  instead  of 
reaction  progress  as  the  second  variable.  The  model  described  here  can  be  used  for  simulations 
of  premixed,  diffusion,  and  partially  premixed  turbulent  combustion  with  high  prediction 
accuracy  and  minimal  computational  cost. 

3.2.1  Model  Description 

The  multi-step  chemistry  model  consists  of  solving  the  fuel  mixture  fraction  variance,  fuel 
reaction  progress  variance,  and  species  transport  equations. 
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The  fuel  mixture  fraction  is  defined  by 


£- 


z-z„. 


zf-zm 


(32) 


where  Zk  is  defined  as  an  atom  that  is  present  in  the  oxidizer  and  not  the  fuel  (or  vice-versa). 
Atomic  nitrogen  is  currently  being  used. 


The  fuel  reaction  progress  is  defined  by 


e = 


v  -Y 

2/.  max  *! 

v  -Y 

x/,  max  x/,min 


(33) 


where  Yr  and  Yt  are  strictly  functions  of  the  fuel  mixture  fraction,  and  Yf  is  the  fuel  mass 

j  ,max  j  ,min  **  j 

fraction. 
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The  variables  £ ,  9,  a ^ ,  and  cr0  are  needed  to  define  the  probability  density  function  (PDF)  that 


is  used  in  determining  PDF-averaged  quantities  of  temperature,  specific  volume,  and  reaction 
rate.  These  averaged  quantities  are  determined  by  integrating  the  instantaneous  values  over  the 


PDF. 


5=  P(0^)d0  d§ 


(34) 


where  a  is  the  quantity  being  averaged  and  p(0,^)  is  the  PDF  as  a  function  of  fuel  mixture 
fraction  and  fuel  reaction  progress,  which  are  assumed  to  be  independent  and  to  be  represented 
as  p{o)p(£).  The  shape  of  the  PDF  is  user-specified.  For  £ ,  the  shape  can  be  a  beta  function  or 
a  top-hat  function.  For  9 ,  the  shape  can  be  one  of  top-hat  or  tri-delta.  The  variances  determine 
the  spread  of  the  PDF  and  are  known  from  solving  the  two  variance  transport  equations  (Eq.  29 
and  30). 


The  instantaneous  reaction  rate  for  step  k  is  obtained  from  the  equation 


=  ATP  Pre 


(-If) 


n- 


(35) 


where  A  is  the  pre-exponential  factor,  /?  is  the  temperature  exponent,  y  is  the  pressure 
exponent,  Ea  is  the  activation  energy,  R  is  the  universal  gas  constant,  c i  is  the  concentration  of 
reactant  j ,  and  a }  is  the  concentration  exponent  for  reactant  j .  The  PDF-averaged  reaction 
rate  for  step  k  is 

7k  =  I\rkp(0)p(g)d0d4  (36) 

The  PDF-averaged  species  net  production  rate  is  determined  from 

k 


where  vi(k)  is  the  stoichiometric  coefficient  for  species  i  in  reaction  k  . 

The  derivative  of  a>i  with  respect  to  Cj  is  needed  for  the  linearized  portion  of  the  species 
transport  equation  source  term.  This  derivative  is  obtained  as 


where 


and 
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(40) 


a :  —  where  c  is  a  reactant  in  step  k 
c  i 

0  where  cy  is  not  a  reactant  in  step  k 

Tire  determination  of  the  functionality  of  the  a(0,%)  term  in  Equation  6  has  not  yet  been  defined. 
Equilibrium,  at  constant  enthalpy  and  pressure,  is  used  to  define  this  functionality.  The 
equilibrium  values  for  temperature,  specific  volume,  and  reaction  rate  all  vary  strongly  as  a 
function  of  fuel  reaction  progress  and  fuel  mixture  fraction.  The  underlying  assumption  in  using 
this  approach  is  that  the  equilibrium  profile  shape  of  a(0,^)  is  the  same  as  the  instantaneous 
profile  shape  of  a  away  from  equilibrium.  Figure  14  shows  a  sketch  of  what  these 
corresponding  profiles  might  look  like  over  fuel  reaction  progress  space. 


Figure  14.  Equilibrium  Concentrations  Versus  Instantaneous  Concentrations  as  a  Function  of 

Reaction  Progress 

A  difficulty  arises  in  the  solution  of  the  transport  equations  when  treating  the  PDF-averaged 
quantities  in  this  manner.  The  averaged  quantity  in  Equation  6,  namely  a ,  is  computed  as  if 
only  a  function  of  fuel  mixture  fraction  and  fuel  reaction  progress.  Fuel  mixture  fraction  and 
fuel  reaction  progress  are  both  functions  of  the  fuel  reaction  step  but  are  not  functions  of  the 
other  steps  in  the  multi-step  global  mechanism.  For  example,  if  the  CO  step  is  far  from 
equilibrium,  then  the  PDF-averaged  quantity  would  be  based  on  the  equilibrium  CO 
concentration  but  the  CO  would  actually  be  far  from  equilibrium.  Thus,  a{0 ,£)  used  in 
Equation  6  would  not  be  dependent  on  the  solved-for  species  concentrations.  A  correction  is 
needed  to  account  for  this  error  in  the  PDF-averaging  approach.  This  correction  consists  of 
using  scaling  parameters  to  correct  for  non-equilibrium  states.  Corrections  need  to  be  made  for 
PDF-averaged  temperature  and  reaction  rates.  Specific  volume  corrections  are  implicitly  made 
due  to  the  temperature  correction.  The  scaling  parameter  consists  of  the  ratio  of  the  solved  for 
value  to  the  equilibrium  value  at  the  mean  values  of  9  and  % . 
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(41) 


This  1 9j  is  then  used  in  correcting  the  reaction  rate  rk  by 

h  -JIM'*  (42) 

j 

The  derivative  rates  are  corrected  using  this  same  ratio  in  this  manner 
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dc 
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The  temperature  is  corrected  in  a  similar  manner 


T’  =  VTe 


eq 


where 


=  1-A  T/ 


eq 


where 


A7’  =  r((fM..JeJ-r((Fy=1...J) 


(43) 


(44) 

(45) 


(46) 


The  two  temperatures  in  Equation  46  are  computed  as  a  function  of  enthalpy  and  PDF-averaged 
species  mass  fractions  and  solved  for  mass  fractions  respectively. 

Having  determined  the  corrected  PDF-averaged  temperature,  specific  volume,  and  species  net 
production  rates,  the  solution  then  focuses  on  solving  the  species  transport  equations.  These 
transport  equations  are  very  sensitive  to  large  source  terms  for  the  consumption  and/or 
production  of  species.  The  species  net  production  rate  is  lagged  one  iteration  to  prevent 
undesired  oscillations  in  the  solution.  Lagging  the  species  source  term  creates  the  drawback  that 
overshoot  might  occur.  This  overshoot  can  cause  species  fractions  to  become  negative  if  not 
somehow  limited  or  relaxed  in  the  solution.  The  current  implementation  limits  the  reaction  rates, 
not  the  species  net  production  rate.  Therefore,  if  negative  species  concentrations  are  being 
solved  for,  then  all  of  the  destruction  steps  for  that  species  are  reduced  by  20%.  This  is 
continued  until  the  species  concentration  is  no  more  than  half  of  its  original  value.  This  limiting 
technique  works  remarkably  well.  Typically,  after  just  a  few  iterations  the  need  for  limiting 
vanishes  for  most  of  the  computational  cells  if  not  for  all  of  them.  In  some  instances  in  the  post 
flame  region,  the  fuel  fraction  is  very  near  zero.  In  these  regions,  limiting  factors  may  still  be 
required  in  order  to  keep  the  solution  well  behaved.  However,  this  type  of  limiting  only  affects 
the  fuel  species  well  downstream  of  the  flame  and  does  not  affect  the  integrity  of  the  overall 
solution. 
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Because  the  PDF-averaged  reaction  rates  are  averaged  over  fuel  mixture  fraction  and  fuel 
reaction  progress  space,  errors  can  arise  in  the  PDF-averaging  near  the  limit  of  <9  =  1.  This  is 
because  of  the  equilibrium  profile  assumed  for  the  PDF-integration.  For  fuel-rich  conditions, 
Y oxidizer eq  ~  0  anc*  f°r  fuel-lean  conditions  Yfueleq  ~  0  .  Thus,  the  fuel  reaction  rate  can  prematurely 


approach  zero.  These  same  errors  occur  near  0  =  0  but  don’t  pose  nearly  the  same  problem  due 
to  the  fact  that  before  the  fuel  has  reacted  the  other  species  in  the  global  mechanism  are  typically 
dormant.  The  PDF-averaged  fuel  reaction  rate  should  never  dip  below  the  laminar  reaction  rate 
in  the  post-flame  region.  Therefore,  in  the  solution  procedure,  when  0  >  0.5  then  the  fuel-step 
reaction  rate  is  treated  as  such: 

=(l-/?)  rPDF  +  /3  7lamn  ar  (47) 


where 


[3  =  max 
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(48) 


This  causes  the  PDF-averaged  fuel  rate  to  slowly  transition  to  the  laminar  rates  as  fuel  reaction 
progress  approaches  one.  Doing  so  greatly  enhances  convergence. 

For  pure  diffusion  flame  problems,  the  reaction  progress  is  sometimes  not  well  defined, 
fact  leads  to  difficulty  in  the  implementation  of  Equation  47.  Instead,  to  overcome  the 
realistically  low  fuel  PDF  rates  the  following  equation  is  used: 

7f  =max(rroF,rtominaf) 

The  species  transport  equations  are  solved  via  an  implicit,  linear  point-solver  method, 
coupled  species  equations  are  solved  at  each  computational  cell.  This  technique  leads  to 
convergence  but  provides  for  a  robust  solution  procedure. 

The  current  species  transport  equation  solution  method  employed  by  the  multi-step  PDF  module 
in  CFD-ACE+  is  inadequate  for  large-scale  simulations  of  combustion  phenomena.  As 
previously  mentioned,  the  existing  method  for  solving  the  species  transport  equations  utilizes  a 
cell-by-cell  point  solver.  Specifically,  the  solver  ‘sweeps’  over  all  of  the  computational  cells  in 
the  domain  solving  the  coupled  species  equations  at  each  cell.  In  a  sense,  this  method  holds  the 
solution  values  at  all  other  cells  constant  while  solving  for  the  species  fractions  at  a  given  cell. 
While  adequate  for  small  problems,  large  problems  require  a  large  number  of  iterations  in  order 
to  assure  global  convergence.  Steady-state  problems  can  usually  afford  a  large  number  of  global 
iterations  to  achieve  ample  convergence.  However,  transient  cases  such  as  LES  need  to  keep  the 
number  of  global  iteration  per  time  step  to  a  minimum  if  practical  run  times  are  to  be  realized. 

The  other  variables  in  the  CFD  calculations  are  solved  for  sequentially  using  a  field  solution 
procedure.  This  method  propagates  boundary  conditions  through  the  solution  domain  very 
efficiently  and  is  the  standard  approach  for  most  CFD  calculations.  Enthalpy  is  solved  for  in  this 
fashion  and  significant  problems  have  been  seen  with  the  coupling  of  the  species  (solved  for 
using  the  point  solver)  and  enthalpy.  Adopting  this  same  strategy  for  the  species  solution  is  a 
plausible  approach  to  achieve  good  coupling  of  the  enthalpy  and  species  and  was  tested,  hi  this 


This 

non- 
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slow 


27 


8503/8 


approach,  a  given  species  was  solved  over  the  entire  computational  domain  while  holding  all  the 
other  species  constant.  The  remaining  species  were  then  solved  for  in  turn  in  a  likewise  manner. 
This  procedure  has  the  benefit  of  fast  propagation  of  boundary  information.  However,  species 
are  highly  coupled  and  exhibit  nonlinear  behavior  and  solutions  using  this  methodology  rarely 
converge. 

The  solution  to  this  problem  would  be  to  modify  the  solver  to  allow  for  a  field  solution  of  the 
coupled  species  equations.  Instead  of  solving  for  species  1  while  holding  species  2,  3,  4,  etc. 
constant,  all  of  the  species  would  be  solved  simultaneously  over  the  entire  field.  In  the  CFD- 
ACE+  solver,  this  would  necessitate  that  the  field  solver  handle  more  than  current  limitation  of 
solving  for  only  one  dependent  variable  at  a  time.  Using  such  a  coupled  field  solution,  better 
coupling  to  other  flow  variables,  such  as  enthalpy,  is  possible.  This  occurs  because  the  species 
information  will  propagate  at  the  same  rate  as  the  enthalpy  information.  Also,  the  species  will 
remain  implicitly  coupled  at  each  cell.  This  solution  procedure  will  implicitly  will  require  more 
CPU  time  per  iteration,  but  will  converge  in  a  fewer  number  of  iterations.  However,  the 
implementing  this  in  CFD-ACE+  was  beyond  the  scope  of  the  present  Phase  II  effort. 

3.2.2  Table  Look-Up 

The  calculation  of  PDF  averaged  species’  chemical  source  terms  with  multi-step  chemistry  is 
extremely  expensive.  A  table  allows  the  integrated  specific  volume,  temperature,  and  reaction 
rates  to  be  stored  as  a  function  of  0 ,  % ,  cre,  cr4,  enthalpy,  and  rj .  To  improve  computational 

efficiency,  this  data  is  first  calculated  and  stored  for  subsequent  table  look-up.  CFDRC  has 
evaluated  several  table  lookup  schemes  for  applicability  with  the  implemented  multi-step  PDF 
combustion  model.  Originally  the  in-situ  tabulation  method  (ISAT),  based  on  the  work  of  Pope 
(1997),  was  proposed  to  represent  compositions  that  are  accessed  during  the  simulation  without 
requiring  storage  for  unneeded  compositions.  For  laminar  flames,  as  chemical  species  evolve 
through  composition  space  along  low-dimensional  manifolds,  table  values  are  needed  for  only  a 
fraction  of  the  allowed  composition  space.  Thus,  a  lot  of  storage  space  is  not  required.  However, 
for  turbulent  combustion  calculations,  the  PDF  averaged  chemical  source  terms  have  a  huge 
memory  requirement.  The  PDF  for  turbulent  combustion  calculations  is  defined  by  the  mixture 
fraction,  mixture  fraction  variance,  reaction  progress  variable,  reaction  progress  variable 
variance,  enthalpy  and  ratio  of  dynamic  to  equilibrium  CO  composition.  Problems  with  ISAT 
were  also  encountered  during  LES  calculations,  due  to  the  resolution  of  the  large-scale 
fluctuations.  The  fluctuations  created  a  large  number  of  potential  points  in  the  ISAT  table  and 
prevented  a  high  level  of  table  retrievals  versus  direct  integrations.  For  the  assumed  PDF 
approach  developed  in  the  SBIR  for  turbulent  combustion  calculations  using  both  RANS  and 
LES,  the  ISAT  approach  did  not  provide  significant  benefits  as  far  as  memory  requirements  and 
computational  speedup  are  concerned.  Instead,  the  relative  expensive  data  search  and 
extrapolation  during  the  data  retrieval  process  made  ISAT  time  consuming. 

In  this  project,  an  alternative  table  look-up  approach  using  an  in-situ  procedure  has  been 
implemented.  A  table  outline  is  created  in  a  pre-processing  stage  with  grid  points  covering  the 
entire  composition  space.  The  interval  spacing  can  be  either  even  or  variable.  Data  at  the  needed 
table  points  is  then  computed  and  stored  as  needed  during  the  simulation.  In  this  way,  only 
chemistry  data  pertinent  to  the  simulation  is  calculated.  Initial  iterations  are  slow,  but  retrieval 
rates  approach  100%  as  convergence  is  achieved.  During  the  simulation  process,  chemical 
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reaction  information  for  each  species  at  each  cell  can  be  obtained  by  checking  the  information 
corresponding  to  the  cell  and  its  neighboring  cells  in  the  table.  This  procedure  is  simpler  than 
the  ISAT  framework  and  achieves  similar  results  in  speeding  computation  times.  LES 
calculations  especially  benefit  from  this  table  approach,  with  retrieval  rates  approaching  100% 
after  a  limit-cycle  condition  has  been  reached. 

3.3  Improved  LES  Inflow  Boundary  Conditions 

One  significant  problem  encountered  in  LES  modeling  of  turbulent  flow  is  the  specification  of 
the  inflow  boundary  conditions.  For  simulations  that  inherently  generate  a  great  deal  of 
turbulence,  specification  of  the  inflow  boundary  condition  turbulence  levels  is  not  that  critical,  as 
the  flow  will  generate  its  own  turbulence.  However,  to  capture  flow  features  that  depend  on  the 
free-stream  turbulence,  such  as  flow  separation  and  boundary  layer  transition,  improved 
boundary  conditions  are  needed.  Research  has  shown  that  application  of  random  turbulence  at 
the  inlet  does  not  capture  the  relevant  length  and  time  scales  of  the  flow  (Schltiter,  2002).  The 
standard  approach  is  to  specify  a  random  inlet  condition  for  each  time  to  give  the  right  values  of 
turbulent  fluctuations  based  on  some  estimated  turbulent  kinetic  energy.  This  is  the  traditional, 
and  simplest  way,  to  handle  specification  of  turbulent  energy  at  the  inlet  boundary  condition. 
However,  the  structures  introduced  with  this  approach  are  small  scale  (on  the  order  of  the  grid 
size)  and  tend  to  dissipate  very  quickly,  with  the  end  result  being  nearly  identical  to  a  laminar 
flow  boundary.  Several  improvements  to  the  LES  modeling  capabilities  with  respect  to  inlet 
boundary  conditions  in  CFD-ACE+  were  completed.  Two  new  approaches  have  been  implement 
to  provide  inlets  with  a  more  representative  turbulent  flow. 

The  first  approach  generates  a  random  inlet  correlated  to  a  specified  length  scale  using  an 
approach  similar  to  that  of  Klein  et  al  (2003).  This  method  gives  the  user  control  of  the  eddy 
size  at  the  boundary  but  does  not  allow  for  a  range  of  sizes.  Therefore,  some  development  length 
is  still  needed  to  allow  the  flow  to  generate  the  full  length  of  length  and  time  scales.  The  second 
method  is  to  simulate  the  inlet  turbulence  via  an  auxiliary  solution  of  the  upstream  region, 
utilizing  periodic  boundary  conditions.  The  resulting  periodic  outflow  boundary  condition 
profile  can  then  be  read  into  the  downstream  simulation  as  an  inflow  boundary  condition.  This 
approach  provides  the  full  range  of  length  and  time  scales  at  the  inlet  boundary,  but  only  works 
in  cases  where  a  simple  channel  flow  can  be  defined  to  characterize  the  inlet  and  is  not  valid  for 
complex  geometries  or  for  inlets  with  swirl. 

The  boundary  conditions  were  tested  using  a  channel  flow  between  two  parallel  plates.  Figure 
15  shows  the  geometry  used  for  the  simulations,  along  with  the  corresponding  grid.  The  top  and 
bottom  boundary  conditions  are  that  of  walls,  while  the  sides  are  simulated  as  periodic. 
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Figure  15.  Geometry  and  Grid  for  Channel  Flow  Between  Two  Parallel  Plates 


If  the  inlet  and  outlet  boundaries  are  taken  to  be  periodic,  then  the  fully  developed  turbulent  flow 
structure  shown  in  Figure  16  results.  The  contours  shown  in  Figure  16  are  of  axial  velocity  in 
the  center  of  the  channel  in  between  the  top  and  bottom  walls.  Figure  17  shows  the  axial 
velocity  at  an  instant  in  time  at  the  periodic  inflow  plane.  The  inlet  velocity  fluctuations  tend  to 
be  distributed  randomly  in  space,  but  there  is  a  definite  range  of  eddy  sizes  in  the  simulations. 
Also,  the  approximate  eddy  length  scales  at  the  inlet  can  be  determined.  The  largest  eddy  is  no 
larger  than  ~  1/3  of  the  channel  height. 
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Figure  16.  Turbulent  Flow  Structures  for  Axial  Velocity  in  Channel  Flow  Between  Two  Parallel 
Plates  Utilizing  Periodic  Inlet-Outlet  Boundary  Conditions 


Figure  1 7.  Inlet  Axial-Velocity  Profiles  for  Periodic  Inlet-Outlet  Simulation  of  Flow  Between 

Two  Parallel  Plates 
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3.3.1  Correlated  Random  Inlet  (CRI) 

Turbulence  has  always  been  characterized  with  randomness  in  the  flow,  so  historically  most  of 
the  Large  Eddy  Simulations  were  performed  by  specifying  a  mean  flow  and  a  random  fluctuation 
over  the  mean.  In  this  approach  (referred  to  as  Gaussian  random),  no  correlation  between 
adjacent  cells  is  assumed.  Although  the  turbulent  flows  have  a  degree  of  randomness  associated 
with  them,  it  has  been  observed  that  the  fluctuations  are  not  completely  random  and  that  a  three- 
dimensional  correlation  exists  between  the  fluctuations.  The  difference  between  the  spatially 
correlated  boundary  and  a  Gaussian  random  boundary  is  that  the  random  number  set  used  to 
characterize  the  velocity  fluctuations  are  related  to  their  neighboring  cell  values  and  also  to  their 
neighboring  temporal  values. 

To  enhance  the  existing  capabilities  of  CFD-ACE+,  three-dimensional  correlations  between  the 
velocity  fluctuations  were  specified  at  the  inlet  boundary  for  Large  Eddy  Simulations.  For 
simplicity  and  low  computational  cost,  a  simple  auto-correlation  function  derived  from  the 
assumption  of  homogenous  turbulence  (Batchelor,  1953)  was  used.  The  shape  of  this  auto¬ 
correlation  function  is: 


K»  =  exp 


(  nr ^ 
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This  functional  dependence  fulfills  some  basic  properties  like  Rull(0)  =  1  and  lim/?uu(r)  =  0 .  The 

r— >co 

Length  scale  ‘L’  for  the  turbulent  flow  has  to  be  specified  at  the  inlet  boundary.  This  scale 
determines  the  eddy  size  observed  at  the  inlet  of  the  computational  domain. 

In  order  to  specify  three-dimensional  velocity  correlations  at  a  two-dimensional  inlet  boundary, 
the  flow  velocities  were  time-correlated  based  on  the  fluctuation  magnitude  and  eddy  size 
specified  at  the  inlet,  as  follows: 


^ eddy 
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Here,  teddy  represents  the  time  for  which  the  fluctuations  are  correlated  at  the  inlet,  Leddy  is  the 
eddy  size  at  the  inlet  and  u  is  the  fluctuation  magnitude  at  the  inlet. 

The  computational  procedure  involved  generating  a  set  of  random  numbers  for  the  whole  inlet 
boundary  faces  and  specifying  the  random  velocity  fluctuation  for  every  computational  face  by 
averaging  over  the  random  numbers  of  adjacent  faces  lying  within  the  range  determined  by  the 
eddy  size.  The  averaging  was  done  by  taking  the  convolution  of  the  auto-correlation  function  and 
the  random  numbers.  As  tire  specified  eddy  size  becomes  small  (less  then  the  grid  cell 
dimension)  the  method  defaults  to  the  Gaussian  random  inlet. 

Profiles  showing  inflow  data  determined  via  a  Gaussian  random  and  three  different  specified 
length  scales  are  shown  in  Figure  18.  For  the  correlated  random  inlet  (CRI),  the  user  specified 
eddy  length  scales  are  0.1,  0.3,  and  2.0  m.  The  channel  height  is  2.0  m  and  the  width  is  1.0m. 
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The  correlation  between  neighboring  cells  needed  to  produce  coherent  fluctuations  at  the  eddy 
size  are  immediately  apparent.  In  the  Gaussian  random  approach,  the  velocity  fluctuations  at 
each  boundary  cell  are  independent  and  the  largest  eddy  sizes  are  very  small.  Equally  important, 
these  fluctuations  randomly  vary  at  each  time  step  and  are  not  correlated  from  one  time  step  to 
another.  This  results  in  very  high  frequency  fluctuations  at  the  inlet,  which  decay  very  rapidly. 


(c)  Leddy  =  0.3  m 


(d)  Leddy  =  2.0  m 


Figure  18.  Inlet  Axial-Velocity  Profiles  Showing  Spatial  Correlation  for  the  Gaussian  Random 

and  SCI  Boundaries 
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As  shown  in  Figure  18,  as  the  user  specified  eddy  size  increases,  the  eddy  size  at  the  inlet 
increases.  Comparing  the  images  to  that  shown  in  Figure  17  for  the  periodic  case,  an  eddy  size 
of  0.3  m  seems  to  most  resemble  the  inlet  profiles.  An  eddy  size  of  0.1  m  yields  too  many  small 
scale  eddies  and  no  large  scale  eddies.  With  an  eddy  size  of  2.0  m,  the  size  of  the  channel 
height,  the  entire  inlet  is  one  eddy  and  the  result  is  similar  to  RANS.  hr  addition  to  being 
spatially  correlated,  the  inlet  profiles  of  Figure  18  are  also  temporally  correlated.  Temporal 
correlation  is  accomplished  through  the  application  of  Taylor’s  hypothesis  to  relate  the  length 
and  time  scales.  The  time  scale  is  equal  to  the  specified  length  scale  divided  by  the  mean 
velocity. 

The  developed  flow  along  the  length  of  the  channel  was  computed  for  the  Gaussian  random  and 
CRI  eddy  sizes  of  0.1  and  0.3  m.  These  results  are  shown  in  Figure  19.  In  the  Gaussian  random 
case  fluctuations  at  the  inlet  quickly  die  out.  The  flow  further  downstream  is  nearly  void  of  any 
turbulent  flow  structure  and  rapidly  laminarizing. 


(a)  Gaussian  Random 


(b)  L  =  0. 1  m 


(c)  L=0.3  m 


Figure  19.  Turbulent  Flow  Structures  for  Axial  Velocity  in  Channel  Flow  Between  Two  Parallel 
Plates  Utilizing  a  Gaussian  Random  Method  for  Generation  of  Inflow  Boundary  Condition  Data 


With  a  specified  eddy  size  of  0.1,  there  is  a  significant  development  region  near  the  inlet  where 
the  small  eddies  dissipate.  However  it  appears  that  larger  eddies  are  forming  downstream  and 
that  with  a  sufficiently  long  development  length,  a  fully  turbulent  flow  profile  would  result.  A 
long  development  length  is  necessary  because  the  specified  size  of  0.1  m  was  too  small  in 
comparison  to  the  fully  developed  flow  modeled  by  the  periodic  inlet-outlet  simulation.  With  a 
specified  inlet  eddy  size  of  0.3  m,  the  development  region  near  the  inlet  is  much  smaller  than  that 
for  an  eddy  size  of  0.1  m.  The  flow  appears  to  be  rapidly  approaching  the  results  shown  in 
Figure  10  for  the  periodic  inlet-outlet  simulation. 
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In  summary,  the  Gaussian  random  inlet  is  very  easy  to  implement  and  has  traditionally  been  used 
as  the  inflow  method  of  choice.  The  main  problem  with  the  Gaussian  random  inlet  is  that  the 
fluctuations  generated  are  not  correlated  in  spatially.  The  eddy  sizes  are  determined  by  the  grid 
and  thus,  the  fluctuations  do  not  represent  true  turbulent  flow.  These  small  scale  fluctuations 
rapidly  dissipate  downstream  of  the  inlet,  basically  resulting  in  laminar  inlet  flow. 

The  correlated  random  inlet  method’s  strength  is  that  the  inlet  eddies  are  not  just  random,  but  are 
correlated  in  both  time  and  space  to  provide  more  realistic  turbulent  profiles.  The  method  is  also 
a  very  inexpensive  calculation,  and  can  be  used  for  swirl  inlets  and  inlets  with  complex 
geometry.  A  disadvantage  of  this  method  is  that  the  user  needs  to  know  the  eddy  length  scales  at 
the  inlet  a  priori.  In  addition,  there  will  exist  a  development  region  downstream  of  the  inlet  as 
the  flow  develops  the  full  range  of  eddy  sizes. 


3.3.2  Auxiliary  Simulation 

The  second  method  for  defining  an  LES  boundary  is  to  read  in  the  boundary  condition  profile  at 
each  time  step  from  an  upstream,  auxiliary  solution  using  periodic  boundary  conditions.  In  the 
generation  of  inflow  data  from  an  auxiliary  upstream  periodic  inlet-outlet  solution,  there  are  two 
elements  that  make  it  possible  to  use  these  results  as  downstream  inflow  boundary  conditions. 
The  first  is  interpolation.  This  allows  the  auxiliary  solution  grid  to  be  different  from  the  grid  that 
the  generated  inflow  boundary  condition  information  is  to  be  applied  to.  In  the  CFD-ACE+ 
algorithm,  the  interpolation  scheme  consists  of  locating  the  four  closest  nodes  and  computing  a 
weighting  factor  based  on  how  far  away  the  nodes  are  from  the  desired  location.  The  weighting 
factor  is  the  inverse  of  the  distance  away  from  the  desired  location,  normalized  by  the  inverse  of 
the  sum  of  all  four  distances.  The  weighting  factors  sum  to  one.  The  value  of  the  dependent 
variable  at  the  desired  location  is  computed  as  the  sum  of  the  product  of  the  weighting  factor  and 
the  dependent  variable  for  the  four  closest  nodes. 


The  second  element  that  makes  it  possible  to  use  an  auxiliary  upstream  periodic  inlet-outlet 
solution  for  a  downstream  inflow  boundary  condition  is  windowing.  Windowing  makes  it 
possible  to  use  a  data  set  consisting  of  a  limited  number  of  frames  (velocity  profiles  at  each  time 
step)  over  and  over  at  the  inlet  without  discontinuities  occurring  when  the  cycle  repeats  itself. 
For  example,  the  auxiliary  solution  might  generate  two  flow  through  times  worth  of  data.  The 
downstream  simulation  then  needs  to  be  able  to  use  these  results  to  simulate  ten  flow  through 
times  in  the  actual  simulation.  The  inflow  data  will  then  be  repeated  five  times.  The  modified 
velocity  at  each  time  step  is  computed  as 


U(t)=U0+(u(t)-U0)sin 


71 


tot  J 


(52) 


where  U(t)  is  the  windowing  modified  velocity,  U0  is  the  velocity  at  the  first  frame  of  the  inflow 
data,  U(t)  is  the  velocity  at  the  nth  frame  of  the  inflow  data,  fn  is  the  number  of  frames  used,  and 
hot  is  the  total  number  of  frames  in  the  data  set.  Using  this  relationship  leads  to  a  smooth 
transition  of  the  velocities  from  the  end  of  the  windowing  frames  to  the  beginning  of  the  next 
repeat  of  inflow  data  while  preserving  the  length  and  time  scales  associated  with  the  original  data 
set. 
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Figure  20  shows  the  channel  flow  simulations  using  inlet  and  outlet  boundaries,  with  the  inlet 
boundary  mapped  from  the  periodic  channel  flow  simulations.  The  velocities  plotted  in  Figure 
20  show  a  very  similar  turbulent  flow  structures  as  those  shown  in  Figure  16.  It  is  interesting  to 
note  that  there  is  no  development  region  for  the  flow  near  the  inlet.  This  is  as  expected,  since  the 
inlet  boundary  condition  information  was  obtained  from  a  similar  simulation  performed  with 
periodic  inlet-outlet  boundary  conditions. 


o  m  o  in  £  g 


Figure  20.  Turbulent  Flow  Structures  for  Axial  Velocity  in  Channel  Flow  Between  Two  Parallel 
Plates  Utilizing  an  Auxiliary  Upstream  Solution  for  Generation  of  Inflow  Boundary  Condition 

Data 


The  approach  of  using  an  auxiliary  upstream  simulation  with  periodic  inflow-outflow  boundary 
conditions  to  generate  turbulent  inlet  profiles  is  a  good  way  to  get  realistic  turbulence  at  the  inlet 
where  the  inlet  geometry  is  such  that  a  channel  flow  can  be  specified.  There  is  no  development 
region  at  the  inlet  because  the  turbulence  provided  is  representative  of  fully  developed  turbulent 
flow.  The  full  range  of  eddy  length  scales  are  determined  as  part  of  the  solution,  thus  requiring 
no  a  priori  determination  of  the  eddy  size.  In  addition,  the  boundary  solution  can  be  interpolated 
onto  different  inlet  meshes  and  still  retain  the  turbulent  flow  structure.  In  theory  it  would  be 
possible  to  use  a  generic  upstream  inflow  boundary  condition  generator  and  then  apply  it  to 
multiple  different  downstream  geometries.  A  disadvantage  of  this  approach  is  that  it  has  the 
potential  to  be  very  expensive.  A  separate  simulation  must  be  preformed  simply  to  get  the 
boundary  condition  information  for  the  real  simulation  of  interest. 

3.4  Global  Mechanism  Generation 

An  accurate  reduced  chemical  kinetic  mechanism  is  needed  to  adequately  predict  flame 
stabilization  and  pollutant  formation  in  advanced  combustors.  Detailed  chemical  kinetic 
mechanisms  for  JP-8-type  aviation  fuels  have  been  under  development  for  some  time  (Dagaut  et 
al.,  1994;  Violi  et  al.,  2002;  Mawid  and  Sekar,  2001)  but  are  limited  in  their  validation  and 
application  at  high  pressure  (20-50  atm).  The  above  cited  mechanisms  by  Dagaut,  Mawid  and 
co-workers  have  been  validated  for  ignition  delay  at  atmospheric  or  low  pressure  conditions. 
Further  work  needs  to  be  done  to  develop  simplified  reaction  mechanisms  that  are  valid  at  high 
pressure. 

The  use  of  these  detailed  chemical  kinetic  mechanisms  for  military  fuels  such  as  JP-8  is  not 
practical  in  three-dimensional  CFD  codes  at  the  present  time.  The  impractical  nature  of  using 
detailed  mechanisms  is  largely  due  to  the  extremely  high  computational  expense  required  to 
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integrate  the  stiff  set  of  chemical  kinetic  equations.  Development  of  computational  procedures 
involving  realistic  kinetic  mechanisms  is  needed,  especially  when  accurate  prediction  of 
emissions  is  needed. 


For  practical  CFD  calculations,  including  those  performed  with  the  multi-step  PDF  model,  a 
reduced  set  of  global  reactions  (~  5-10  steps  for  JP-8,  including  NOx)  that  use  simple  Arrhenius 
expressions  for  each  reaction  rate  is  needed.  In  this  SBIR  program,  CFDRC  has  developed  a 
global  multi-step  mechanism  for  JP-8  applicable  to  gas  turbine  combustion  simulations. 
Examples  of  the  global  mechanisms  similar  to  what  has  been  developed  for  JP-8  include  the 
reaction  mechanisms  of  Nichol  et  al.  (1997),  Dryer  and  Glassman  (1973),  Westbrook  and  Dryer 
(1984),  and  Jones  and  Lindstedt  (1988).  The  Nichol  3-step  mechanism  for  methane-air 
combustion  is  shown  in  Table  4. 


Table  4.  The  Nichol  (1995)  Mechanism  for  CH4 


CH,  +-0^C0  +  2H20 
2  - 

co+^o2=>co2 

co2^>co+^o2 


R  =  1011'61  exp 


— 15866^1 


R  =  IO1400  exp 


R  =  IO469  exp 


v  T  ) 

r-i68n' 

V  T  J 

f-  20773 ^ 
T  , 


I \[cH4r[o2r 

[corio2r 

[co2r 


The  global  mechanisms  are  very  easy  to  use  and  cost  effective,  due  to  the  simple  Arrhenius  rate 
expressions  (see  equation  25).  However,  the  mechanisms  are  often  only  applicable  at  or  near  the 
conditions  where  the  experimental  data  used  in  determining  the  Arrhenius  rates  were  collected  In 
order  to  more  accurately  predict  species  and  temperature  profiles  in  practical  CFD  simulations, 
global  multi-step  mechanisms  with  Arrhenius  rates  have  been  developed  in  this  Phase  II  effort. 
CFDRC  has  implemented  a  computational  tool  to  automatically  generate  the  multi-step  global 
mechanism  and  Arrhenius  rates  specifically  for  each  simulated  flow  condition. 

r  -  A- TnPme **  [C, ]0' \C2 ]fl2  (53) 

Global  mechanisms  using  standard  Arrhenius  rates  have  several  significant  advantages.  With 
Arrhenius  rates,  the  source  terms  can  be  directly  evaluated  and  the  solution  of  a  separate  set  of 
equations  is  not  required.  Also,  linearization  of  reaction  rate  source  terms  for  CFD  codes  is  very 
simple  with  Arrhenius  rates  (shown  in  Equation  4),  since  an  analytic  derivation  of  the  reaction 
rate  is  easily  formed.  In  multi-step  mechanisms,  steps  include  the  fuel  oxidation  step  to  CO, 
along  with  one  more  CO  oxidation  step.  Other  steps  accounting  for  NOx,  or  other  minor  species, 
can  be  added  if  needed.  Also,  if  unbumed  hydrocarbons  or  smoke  calculations  are  required,  fuel 
intermediates  may  be  included.  The  coefficients  in  the  Arrhenius  rate  expressions  are 
determined  through  curve  fitting  to  experimental  data. 
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3.4.1  Approach 

The  approach  that  CFDRC  is  following  in  this  SBIR  program  utilizes  zero-  and  one-dimensional 
flow  calculations  with  a  detailed  kinetic  mechanism.  Examples  are  Perfectly  Stirred  Reactor 
(PSR),  Plug  Flow  Reactor  (PFR),  and  laminar  one-dimensional  premixed  flame  calculations. 
Detailed  kinetics  are  used  in  these  models  to  determine  temperature  and  species  profiles  as  a 
function  of  either  residence  time  or  distance.  An  optimization  routine  (Steepest  Descent  method) 
is  then  used  to  fit  global  kinetic  parameters  that  will  match  these  profiles.  Use  of  these  approach 
to  determine  global  kinetic  parameters  is  commonly  referred  to  as  Chemical  Reactor  Modeling 
(CRM).  Figure  21  shows  a  flow  chart  of  this  process.  Under  this  SBIR,  CFDRC  has  developed 
the  capability  to  efficiently  generate  a  reduced  global  mechanism  with  optimum  Arrhenius  rates, 
similar  to  that  shown  in  Table  7,  for  a  given  set  of  conditions.  The  software  is  designed  to 
automatically  determine  the  optimum  coefficients  in  the  Arrhenius  reaction  rates  based  on  CRM 
calculations  using  the  full  detailed  mechanism.  Nicol  and  Make  (1996)  successfully  used  this 
type  of  approach  on  a  limited  basis 


Finished  Results 
to  Flow  Solver 


Figure  21.  Flow  Chart  of  the  Process  Used  to  Automatically  Generate  Reduced  Global 

Mechanisms  from  Detailed  Mechanisms 


3.4.2  PSR-PFR  Combinations 

A  gas  turbine  combustor  can  be  represented  by  various  combinations  of  perfectly-stirred  reactors 
(PSR)  and  plug-flow  reactors  (PFR).  This  type  of  approach  is  frequently  used  to  model  complex 
combustion  systems  (Turns,  2000)  and  has  been  coined  “chemical  reactor  modeling”,  or  CRM. 
One  possible  configuration  is  a  PSR  followed  by  a  PFR,  as  was  used  by  Nicol  (1995).  Figure  22 
shows  this  configuration. 


The  PSR  is  well  mixed,  and  has  not  yet  approached  equilibrium.  This  stage  represents  the 
primary  zone  of  the  combustor,  where  good  recirculation  enhances  the  mixed-ness.  The  PFR  is 
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representative  of  the  secondary  and  dilution  zones  where  mixing  is  not  as  strong  and  the  kinetics 
approach  equilibrium.  In  the  PSR,  the  fuel  reacts  mainly  to  CO;  while  in  the  PFR  section,  the 
CO  is  slowly  oxidized  to  CO2.  Although  this  type  of  PSR-PFR  combination  is  a  simplistic 
representation  of  the  complexities  of  gas  turbine  combustion,  it  is  very  amenable  to  detailed 
kinetic  modeling.  Multidimensional,  high  fidelity  CFD  models,  while  accurately  capturing  the 
relevant  flow  patterns,  are  not  equipped  to  handle  the  coupling  of  detailed  kinetic  calculations. 

One  use  for  CRM  is  to  generate  rate  parameters  for  simple  multi-step  global  mechanisms  with 
Arrhenius  rates  that  can  then  be  used  in  CFD  models  for  gas  turbines.  The  first  step  in  the 
approach  for  obtaining  global  Arrhenius  rate  parameters  is  to  perform  CRM  using  a  detailed 
mechanism.  The  detailed  calculations  were  performed  in  the  PSR  and  PFR  combination,  and  the 
species  net  production  rates  are  generated  as  a  function  of  reactor  residence  time.  The  generation 
of  these  production  rates  can  be  performed  for  single  or  multiple  equivalence  ratios,  at  a  given 
inlet  temperature,  and  for  a  given  operating  pressure. 

Cantera,  an  open  source  code  for  modeling  chemical  reactions  developed  by  David  Goodwin 
(2004),  was  used  as  the  method  of  choice  for  obtaining  the  PSR-PFR  results.  The  PSR 
calculations  were  performed  at  a  residence  time  just  above  the  blow-out  residence  time.  This 
was  determined  in  Cantera  by  systematically  reducing  the  residence  time  to  the  point  where  if 
lowered  any  further,  reactions  cease  to  occur.  This  is  known  as  the  blow-out  residence  time. 
The  blow-out  residence  time  is  used  for  the  PSR  for  two  reasons.  First,  running  the  PSR  just 
above  the  blow-out  condition  allows  for  the  largest  possible  temperature  range  and  thus  makes 
the  global  mechanism  valid  over  this  temperature  range.  Secondly,  at  the  blow-out  residence 
time  the  CO  fraction  is  near  its  peak  value,  thus  providing  a  good  range  of  CO  concentration 
over  which  to  fit  the  global  mechanism. 

The  PFR  calculations  were  initiated  using  the  exit  conditions  of  the  PSR  calculation.  The  PFR 
was  then  allowed  to  react  over  a  range  of  residence  times  that  approached  equilibrium.  Results, 
consisting  of  temperature,  species  net  production  rates,  and  species  concentrations,  are  recorded 
over  the  range  of  residence  times  and  again  over  the  range  of  equivalence  ratios  involved. 

Next,  the  temperature  and  species  concentrations  at  each  residence  time  are  used  with  the  global 
mechanism  to  determine  the  net  production  rates  predicted  by  the  global  mechanism. 
Throughout  the  PSR-PFR  combination,  temperature  and  species  concentrations  are  varying,  thus 
providing  a  wide  range  of  parameters  over  which  the  global  net  production  rates  are  calculated. 

The  goal  is  then  to  match  the  globally  computed  production  rates  with  the  production  rates 
predicted  by  the  detailed  mechanism.  Therefore,  the  difference  between  the  global  rates  and  the 
detailed  rates  needs  to  be  minimized.  Standard  optimization  techniques  (namely  the  steepest 
descent  method)  are  used  to  minimize  the  error  by  varying  the  global  Arrhenius  rate  parameters. 
This  optimization  routine  has  been  embedded  in  the  routines  for  computing  the  PSR-PFR  results. 

Although  the  rates  may  match  very  well,  what  is  of  real  interest  is  the  match  between  the  species 
concentrations  and  temperature  computed  from  the  global  mechanism  for  the  same  conditions 
used  for  the  detailed  mechanism  calculations.  To  demonstrate  applicability,  the  global 
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mechanism  is  then  used  to  compute  temperature  and  species  profiles  for  the  same  PSR-PFR 
configuration. 

In  order  to  facilitate  the  computations  of  the  global  mechanism  based  PSR-PFR  combination,  a 
separate  code  was  required  that  could  deal  with  global  mechanisms.  This  code  was  written  using 
the  PETSc  solvers  (Balay  et  al.,  2001).  The  PSR  equations  were  solved  using  a  Newton- 
Raphson  non-linear  solve  with  pseudo  time  stepping.  The  PFR  equations  were  solved  also  using 
a  Newton-Raphson  non-linear  solve  using  the  implicit  backwards  Euler  method  for  time 
stepping.  Cantera  was  used  for  the  property  calculations:  enthalpy,  density,  mean  molecular 
weight,  etc. 

Some  examples  of  the  types  of  mechanisms  that  can  be  generated  are  shown  below,  both  for 
methane  and  JP-8.  These  examples  aim  to  show  the  flexibility  of  the  CRM  approach  and  are  not 
intended  to  be  the  only  or  even  the  best  multi-step  global  mechanisms  possible  for  these 
conditions.  Using  the  newly  developed  CRM  capabilities,  mechanisms  can  easily  be  generated 
for  any  conditions  of  interest  for  both  methane  and  kerosene  fuels. 

Methane  3-Step 

A  3-step  global  mechanism  was  generated  for  methane  oxidation.  The  GRI-3.0  mechanism, 
consisting  of  53  species  and  325  reactions,  was  used  for  the  detailed  kinetic  calculations.  The  3- 
step  global  mechanism  is  of  the  same  form  as  that  of  Westbrook  and  Dryer  (1984).  The 
optimized  Arrhenius  parameters  are  shown  in  Table  5.  These  parameters  are  valid  over  the 
equivalence  ratio  range  of  0.6- 1.0,  an  inlet  temperature  of  600  K,  and  an  operating  pressure  of  1 
atm. 


Table  5.  Optimized  Arrhenius  Rate  Parameters  for  a  Global  3-step 
Methane  Oxidation  Mechanism 


CH4  +~02  =>  C0  +  2H20 

R{  =7.2x10 10  r0177  exp 

^  — 19172^ 

v  T  J 

[ch„]  “*2[oJ  0,“° 

co  +  ^o2  ^>co2 

R2  =3.9x10 n  Tom  exp 

f-9771Vonn[o,r™ 

\  T  ) 

co2=>co+^o2 

=1.8x106T_0121  exp 

(- 15852^ 

l  T  j 

[CO,]"192 

The  comparison  of  the  detailed  net  production  and  global  net  production  at  the  equivalence  ratio 
of  0.6  is  shown  below,  Figure  23,  as  a  function  of  reactor  residence  time.  The  initial  point  is 
determined  from  the  PSR  calculation.  The  remaining  points  are  from  the  PFR  simulation.  The 
global  mechanism  does  a  remarkably  good  job  in  approximating  the  net  production  rates  for  the 
species  involved. 
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Figure  23.  Comparisons  of  Predictions  in  Net  Species  Productions  for  Detailed 

and  Global  Mechanisms 

Additionally,  the  temperature  and  species  mole  fractions  predicted  by  the  global  mechanism  are 
compared  to  those  calculated  via  the  detailed  mechanism.  Figures  24  and  25  show  these 
comparisons.  It  can  be  seen  that  initially,  the  temperature  is  slightly  higher  for  the  global 
mechanism  at  the  exit  of  the  PSR.  However,  the  temperature  quickly  approaches  that  of  the 
detailed  mechanism.  The  species  concentrations  are  captured  fairly  accurately  except  for  a  slight 
offset,  which  is  most  likely  due  to  the  limited  amount  of  species  contained  in  the  global 
mechanism.  Water,  in  particular,  seems  to  deviate  from  the  detailed  calculation.  Most  likely  a 
global  mechanism  involving  more  steps  would  be  needed  to  accurately  predict  the  water  mole 
fraction. 


Figure  24.  Comparison  of  Temperature  and  Species  Profiles  Computed  via  Global  and  Detailed 
Mechanisms  for  <jr=0.6,  Tin=600  K,  P=1  atm,  and  a  PSR  Residence  Time  of  0.1  msec 
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PFR  Residence  Time  (sec)  P1R  Residence  Time  (sec) 

Figure  25.  Comparison  of  Species  Profiles  Computed  via  Global  and  Detailed  Mechanisms  for 
<p=0.6,  Ti„=600K,  P=1  atm,  and  a  PSR  Residence  Time  of  0.1  msec 


.TP-8  5-Step 

A  detailed  mechanism  for  JP-8  combustion  which  consisted  of  200  species  and  1345  reactions 
was  used  for  the  detailed  kinetic  calculations.  From  this  detailed  mechanism,  two  5-step  global 
mechanisms  were  generated  for  JP-8  combustion  (Tables  6  and  7).  Both  involve  an  intermediary 
step  containing  H2.  The  first  mechanism  (Table  6)  utilizes  the  water-gas  shift  reaction.  The 
second  mechanism  (Table  7)  involves  C12H23  being  oxidized  to  CO  and  H2.  CO  and  H2  then  are 
further  oxidized  by  O2  in  steps  2  and  3.  The  rate  parameters  for  both  mechanisms  are  shown 
below  in  the  table. 


Table  6.  Global  5-Step  Mechanism  for  C12H23  Combustion  Using  a  Water-Gas  Shift  Reaction 


C,  1H73+602  =>  12CO  +  11. 5// 2 

R,  =  3.24  xlOI2r0100  exp 

r'14r798)c«<]l0,J[oJ“,! 

CO  +  H20  =>  C02  +  H2 

R2  =3.3xl06ra'°09exp' 

\ 

-4242\ 

T  J 

CO] 0  666  [o2] 0730 

co2  +  h2^co  +  h2o 

/ 

R3  =  7.01  xl07r° 148  exp 

\ 

— 1564 1^1 

.  T  ) 

[co2  ]  0  774  [//2] 0542 

h2  +  -^o2  => h2o 

R4  =7.08xl012r-0039  exj 

(  -  4446 
\  T 

\h2]'«2[o2Ym' 

J 

h2o^h2  +io2 

R5  =  1.57  xl067”a132  exp 

(-1«70)[h2o]'"’ 
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Table  7.  Global  5-Step  Mechanism  for  C12H23  Combustion  Using  a  CO  Oxidation  Step 


Cl2H23  +60,  =>12C0  +  11.5//2 

R,  =  6.44xlO“r0075  exp^ 

\ 

— 13593^ 
.  T  ) 

[c//4]  0  94,[o2]  0930 

co+^o2=>co2 

R2  =3.76 xlOuT~om  exp 

f-7955> 

* 

[co]1  392  [o2f 240 

T 

V  1  J 

co2=>co+^o2 

R3  =  1.72  xl06r'0136  exp 

"-15236" 

[C02]'ni 

V  T  , 

h2  +  ^o2=>h2o 

Ra  =  5.75  x  1013  T-0'065  exp 

'150  V  1 1.903  fy-v  11.505 

|  T  yi2\  [o2\ 

h2o^>h2  +io2 

R5  =5.00x1 057_0  172  exp 

15564' 

[h2o]U54 

v  T  , 

Both  mechanisms  were  curve  fit  over  the  equivalence  ratio  range  of  0.6- 1.0,  a  pressure  of  1  atm, 
and  an  inlet  temperature  of  600  K.  Comparisons  between  the  net  species  production  rates 
predicted  by  the  two  mechanisms  at  an  equivalence  ratio  of  0.6  are  shown  in  Figures  26  and  27. 
From  these  figures,  it  can  be  seen  that  the  second  mechanism,  the  one  not  involving  the  water- 
gas  shift  reaction,  seems  to  better  approximate  the  rates  computed  from  the  detailed  mechanisms. 
However,  the  first  mechanism  approximates  the  formation  of  water  more  accurately.  Overall, 
good  agreement  is  shown  between  the  detailed  and  global  mechanisms. 


Figure  26.  Species  Net  Production  Comparisons  for  Detailed  and  Global  Mechanism 

(Mechanism  1)  Calculations 
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Figure  27.  Species  Net  Production  Comparisons  for  Detailed  and  Global  Mechanism 

(Mechanism  2)  Calculations 

One  drawback  of  this  approach  for  determining  the  Arrhenius  parameters  is  that  matching  the 
species  net  production  rates  does  not  necessarily  mean  that  the  temperature  profile,  nor  the 
species  concentrations  will  match.  An  improvement  on  this  algorithm  would  be  to  minimize  the 
difference  between  the  temperature  and  species  concentrations  computed  from  the  detailed  and 
global  mechanisms.  Thus,  the  net  production  rates  would  not  match,  but  the  quantities  of  interest 
would.  This  type  of  optimization  is  much  more  involved  however,  as  it  requires  performing 
PSR-PFR  calculations  with  the  global  mechanism  at  each  iteration  of  the  optimization  solver. 

3.4.3  Laminar  One-Dimensional  Flame  Calculations 

Another  approach  to  obtaining  global  rate  parameters  is  to  use  detailed,  laminar,  one¬ 
dimensional,  freely-propagating  flame  calculations.  This  approach  avoids  several  drawbacks  of 
the  PSR-PFR  approach.  Namely,  the  PSR-PFR  approach  ‘ignores’  the  temperature  range  from 
the  inlet  temperature  to  the  PSR  temperature.  Typically,  this  is  a  temperature  range  in  excess  of 
1000  K.  Secondly,  the  PSR-PFR  approach  doesn’t  guarantee  that  the  laminar  flame  speed  will 
be  accurately  predicted.  The  laminar  1-D  calculations  cover  the  entire  temperature  range  (i.e. 
inlet  to  equilibrium).  Also,  the  laminar  flame  speed  is  inherent  in  the  calculation. 

There  are  several  steps  involved  in  obtaining  global  rates  from  detailed  1-D  laminar  calculations. 
These  are  (1)  solve  the  1-D  freely  propagating  flame  using  the  detailed  kinetic  mechanism,  (2) 
select  and  assemble  data  that  is  to  be  fit  against,  and  (3)  perform  optimization  of  Arrhenius 
parameters  using  the  global  mechanism  in  the  1-D  calculation.  Each  of  these  steps  will  be 
discussed  in  detail. 

First,  the  1-D  freely  propagating  flame  is  solved  with  Cantera.  The  detailed  mechanism  is  used 
in  this  calculation.  This  involves  an  iterative  process  of  finding  the  inlet  velocity  at  which  the 
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flame  detaches  from  the  inlet  and  moves  downstream  and  sets  up  a  stable  flame  with  zero 
temperature  gradients  on  either  side. 

Second,  the  data  (really  the  detailed  mechanism  results)  desired  to  fit  the  global  mechanism 
against  needs  to  be  assembled.  This  can  be  temperature,  species  (those  that  exist  in  the  global 
mechanism)  fractions,  or  species  net  production  rates.  Any  combination  of  these  can  be  used. 
For  the  PSR-PFR  model  discussed  earlier,  only  species  net  production  rates  were  used.  For  the 
current  laminar  1-D  calculations,  temperature  and  species  are  used.  However,  any  combination 
of  temperature,  species  fractions,  and/or  species  net  production  rates  could  be  used. 

Finally,  the  optimization  is  performed  using  initial  guesses  for  the  global  rate  parameters.  A  1-D 
flame  is  computed  via  Cantera  using  these  global  kinetics.  The  ‘wellness’  of  the  fit  is  then 
calculated.  For  example,  this  could  be  the  difference  between  the  temperature  profile  calculated 
by  the  detailed  mechanism  and  the  global  mechanism.  The  derivatives  of  this  ‘function’  with 
respect  to  the  Arrhenius  parameters  are  then  computed.  This  involves  multiple  1-D  flame 
calculations  and  can  become  very  expensive  for  a  large  number  of  iterations.  At  each  iteration,  a 
new  set  of  parameters  that  give  a  lower  value  of  the  fitting  function  are  found.  This  procedure 
continues  until  the  local  minimum  is  reached,  and  at  this  point  the  optimum  Arrhenius  kinetic 
parameters  for  the  global  mechanism  are  found. 

Methane  3-Step 

To  demonstrate  this  approach,  a  three-step  global  mechanism  was  developed  for  premixed 
methane  combustion.  The  GRI-Mech  3.0  was  used  for  the  detailed  mechanism.  This  mechanism 
consists  of  53  species  and  325  reactions.  The  global  mechanism  consisted  of  5  species  and  3 
reaction  steps.  This  mechanism  was  fit  at  the  operating  conditions  of  <(>=1.0,  T=300  K,  and 
P=1 .0  atm.  The  curve  fit  parameters  are  shown  below. 


Table  8.  Global  3-Step  Mechanism  for  Methane  Combustion  at  <(>=1.0,  T=300  K,  and  P=1  atm 


2 CH4  +  30,  =>  ICO  +  4 H20 

R,  =  4.57  xl09r“°'"9  exp 

f -13226 ^ 

l  T  , 

[cha]°***[o2]  0  9221 

2  CO  +  02  =>  2  C02 

Rz  =  3.73xl0117’t)'166  exp 

^  — 12127 

[cg]1843[o2]09345 

v  T  J 

2  C02  =>  2  CO  +  02 

R3  =  6.12x1057’“°'352  exp 

(-15519" 

l  T  , 

[cop1 

Plots  of  temperature  and  species  comparisons  are  shown  in  Figures  28a  and  28b,  respectively. 
The  global  mechanism  optimized  via  the  laminar  1-D  calculation  matched  the  laminar  flame 
speed  almost  exactly.  The  fitness  of  the  optimization  for  this  mechanism  was  based  on  the 
temperature  profile  only.  Therefore,  as  shown  in  Figure  28a,  the  temperature  profile  fit  is 
extremely  good.  As  for  the  species  profiles,  the  general  trend  is  captured,  but  the  dynamics  of 
the  species  profiles  are  not  captured  very  well. 
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(a)  (b) 

Figure  28.  (a)  Temperature  and  (b)  Species  Plots  Showing  Comparisons  Between  Detailed 

and  Global  Calculations 

In  order  to  obtain  a  more  accurate  fit  of  the  species  profiles,  another  mechanism  was  generated. 
This  time,  the  species  profiles  were  included  in  the  optimization.  Weighting  factors  were  applied 
to  each  fitted  profile.  Temperature  had  a  weight  of  1,  CO  had  a  weight  of  0.2,  and  CO2  had  a 
weight  of  0.02.  The  resulting  mechanism  is  shown  below  in  Table  9.  It  is  worth  noting  that  this 
mechanism  and  the  previous  mechanism  only  differ  slightly.  This  illustrates  how  sensitive  the 
solution  is  to  the  kinetic  parameters. 


Table  9.  Global  3-Step  Mechanism  for  Methane  Combustion  at  <|>=1.0,  T=300  K,  and  P=1  atm 


2CH4  +  3 02  =$  ICO  +  4 H20 

R{  =4.56  xl09r"° 0711  exp 

f-14014' 

’  T 

V  1  2 

[c//4]05819[o2]  09223 

2  CO  +  02^>  2  C02 

R2  =3.78xl011r°'221  exp' 

\ 

~923l)[co]L*22[o2]  0  9280 

,  T  J 

2C02  =>  2CO  +  02 

R3  =6.07  xl05r'0108  exp 

f  —  1 5702 

l  T  ) 

[co2]'01i 

Comparisons  of  temperature  and  species  profiles  are  shown  below  in  Figures  29a  and  29b, 
respectively.  The  improvement  in  the  species  fit  is  dramatic.  Now,  the  species  profiles  not  only 
match  the  general  trend,  but  also  match  the  peaks  and  curvature.  The  sacrifice  for  this  match  is 
shown  in  the  temperature  profile.  Whereas  before  the  temperature  profile  matched  almost 
perfectly,  now  there  is  a  slight  deviation  downstream  of  the  initial  flame-front.  Weighting  allows 
the  user  to  focus  the  profiles  on  the  ones  of  most  interest. 
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(a)  (b) 

Figure  29.  (a)  Temperature  and  (b)  Species  Plots  Showing  Comparisons  Between  DetaileD 

and  Global  Calculations 

In  addition  to  the  mechanisms  described  above,  any  extension  can  be  made.  For  example,  if  a 
global  NOx  step  is  desired,  it  can  be  included  and  the  mechanism  can  be  re-optimized  while 
applying  the  appropriate  weighting  factor  for  the  NOx  profile.  If  intermediates,  like  OH  or  H,  are 
desired  then  steps  for  these  species  can  be  included  as  well. 

3.4.4  Opposed  Diffusion  Flame  Calculations 

The  capability  to  compute  global  rate  parameters  from  1-D,  laminar,  opposed  diffusion  flame 
calculations  has  been  implemented.  This  technique  is  similar  to  that  for  the  1-D  laminar 
adiabatic  flame.  However,  in  this  flame  fuel  enters  from  one  side  and  oxidizer  enters  from  the 
other.  Figure  30  shows  the  general  configuration  of  an  opposed  diffusion  flame.  The  drawback 
of  the  freely  propagating,  laminar,  adiabatic  flame  is  that  only  one  equivalence  ratio  is  taken  into 
account.  This  drawback  is  eliminated  in  the  opposed  diffusion  flame.  In  this  flame,  the 
equivalence  ratio  goes  from  zero  at  the  oxidizer  inlet  all  the  way  to  infinity  at  the  fuel  inlet. 
Thus,  the  entire  range  from  fuel  rich  to  fuel  lean  is  covered. 


Opposed  Diffusion  Flame 
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Using  opposed  diffusion  flame  calculations,  rate  parameters  for  a  3-step,  6-species  global 
mechanism  were  generated.  Initially,  the  mechanism  was  fit  to  the  species  net  production  rates. 
The  results  are  listed  in  Table  10.  Finally,  these  rates  were  fit  to  the  temperature  and  species 
profiles  for  a  1-D  laminar  opposed  diffusion  flame  (Table  1 1). 

Table  10.  Intermediate  Arrhenius  Parameters  Optimized  On  Species  Net  Production  Rates 


Reaction 

A 

P 

Ea 

(J/kmol) 

Reactant 

Order 

Reactant 

Order 

2  CR)  +  3  02  => 

2  CO  +  4  H20 

5.46E+10 

0.413 

1 .80E+08 

CH4 

0.4080 

o2 

1.305 

2  CO  +  02  =>  2  C02 

9.96E+06 

-0.203 

1.11  E+08 

02 

0.5356 

CO 

0.7947 

2  C02  =>  2  CO  +  O2 

1.84E+06 

-0.121 

1.78E+08 

co2 

0.5813 

Table  1 1.  Arrhenius  Parameters  Optimized  On  Temperature  And  Species  Profiles 


Reaction 

A 

p 

Ea 

(J/kmol) 

Reactant 

Order 

Reactant 

Order 

2  CH4  +  3  02  => 

2  CO  +  4  H20 

3.51E+10 

0.515 

1.77E+08 

CH4 

0.5662 

02 

1.263 

2  CO  +  O2  =>  2  CO2 

8.24E+06 

-0.005 

1.01  E+08 

o2 

0.5338 

CO 

0.657 

2  C02  =>  2  CO  +  02 

3.19E+06 

-0.346 

1.78E+08 

co2 

0.5460 

The  intermediate  rate  parameters  (Table  4)  are  remarkably  close  to  the  final  ones  (Table  5).  This 
indicates  that  the  rate  parameters  obtained  from  the  species  net  production  rates  are  a  very  good 
starting  point  for  fitting  to  temperature  and  species  profiles.  The  determination  of  the  rate 
parameters  in  Table  is  a  very  inexpensive  calculation  (1-2  orders  of  magnitude  less),  whereas 
fitting  to  temperature  and  species  profiles  is  a  more  costly  calculation.  Thus,  by  obtaining  the 
rates  in  this  two-step  procedure  the  overall  cost  of  computing  the  rate  parameters  is  reduced. 

Figures  31  to  33  show  the  results  of  the  fitted  global  mechanism  compared  to  the  detailed 
mechanism  predictions.  Figure  31  shows  the  temperature  comparisons,  while  Figures  32  and  33 
compare  species  profiles.  Other  than  a  slight  overshoot  in  the  peak  temperature,  the  match  is 
very  good.  The  species,  especially  CO  and  CO2,  match  extremely  well  for  such  simplified 
chemistry. 
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Figure  31.  Temperature  Profile  For  1-D  Opposing  Diffusion  Flame  Using  Fitted  Global  Rate 
Parameters  Compared  To  The  Detailed  GRI-3.0  Mechanism  Temperature  Profile 


Distance  (m) 


O  CH4  —  Detailed 

-  CH4  -  Global 

A  02  —  Detailed 
02  -  Global 
+  N2  —  Detailed 
-  N2  -  Global 


Figure  32.  Species  Profiles  Of  CH 4,  Oj,  And  N2  For  1-D  Opposing  Diffusion  Flame  Using 
Fitted  Global  Rate  Parameters  Compared  To  The  Detailed  GRI-3.0  Mechanism  Species  Profiles 
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Figure  33.  Species  Profiles  Of  CO,  H2O,  And  CO2  For  1-D  Opposing  Diffusion  Flame  Using 
Fitted  Global  Rate  Parameters  Compared  To  The  Detailed  GRI-3.0  Mechanism  Species  Profiles 

Two  additional  steps  were  added  to  the  global  mechanism  previously  described  to  account  for 
NOx  formation.  This  resulting  mechanism  is  shown  in  Table  12.  In  the  first  NOx  step,  the  rate  is 
not  only  a  function  of  the  N2  and  O2  concentrations,  but  also  of  CO2.  In  the  second  step,  the  rate 
is  a  function  of  only  O2  and  N2.  The  first  step  accounts  for  non-thermally  formed  NQ  via 
prompt,  nitrous  oxide,  NNH,  and  non-equilibrium  Zeldovich  pathways.  This  is  because  CO  is 
indicative  of  the  free-radical  chemistry  that  is  linked  to  non-thermal  NO  formation  (Nicol  et  al., 
1997).  Thermal  NOx  is  accounted  for  by  the  second  reaction.  Figure  34  shows  the  comparison 
of  the  NO  detailed  profile  with  that  of  the  global  calculation.  The  sharp  nonlinear  behavior 
around  the  peak  is  not  correctly  captured,  but  the  overall  trend  is  predicted  quite  well.  The  other 
species  matches  are  very  similar  to  those  previously  shown  in  Figure  32. 


Table  12.  Global  Mechanism  Containing  5  Steps  And  7  Species  That  Includes  NOx  Chemistry 


Reaction 

A 

b 

Ea 

(J/kmol) 

Order 

Order 

Order 

2  CH4  +  3  o2  => 

2  CO  +  4  H20 

3.46E+10 

0.515 

1.77E+08 

CH4 

0.4974 

o2 

1.284 

2  CO  +  1  02  => 

2  C02 

1.01E+07 

0.0141 

1.01E+08 

02 

0.5104 

CO 

0.6528 

2  C02  => 

2  CO  +  1  02 

3.87E+06 

-0.336 

1.77E+08 

C02 

0.5316 

1  n2  +  1  o2  +  1  CO  => 

2  NO  +  1  CO 

1.75E+09 

-0.132 

1.08E+08 

02 

1.144 

CO 

1.176 

n2 

0.2348 

1  N2+  1  O2  => 

2  NO 

1.15E+08 

-0.334 

1.21E408 

02 

1.35 

n2 

1.236 
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Figure  34.  Species  Profiles  Of  CO,  H20,  and  CO 2  for  1-D  Opposing  Diffusion  Flame  Using 
Fitted  Global  Rate  Parameters  Compared  to  the  Detailed  GRI-3.0  Mechanism  Species  Profiles 


3.5  Initial  Code  Validation 

Comprehensive  validation  is  a  critical  aspect  of  the  program  to  demonstrate  the  capabilities  of 
the  multi-step  chemistry  models  being  developed  in  this  SB1R.  Initially,  two  validation  cases 
have  been  identified  to  test  the  performance  of  the  new  CFD  capabilities  developed  in  this 
program.  A  turbulent  premixed  and  a  diffusion  flame  was  selected  to  demonstrate  the  wide 
applicability  of  the  approach  to  the  analysis  of  turbulent  reacting  flows.  The  following  sections 
discuss  the  two  cases  and  the  CFD  simulations  performed  in  the  SBIR  program. 

3.5.1  University  of  Sydney  Diffusion  Flame  Diffusion  Flame 

The  multi-step  PDF  model  was  first  validated  using  a  diffusion  flame.  The  piloted  burner  tested 
at  the  University  of  Sydney  (Masri  and  Pope,  1990)  has  a  central  jet  of  fuel,  7.2  mm  in  diameter 
surrounded  by  a  stoichiometric  premixed  and  fully  burned  annulus  of  C2H2,  H2,  and  air  with  the 
C/H  ratio  adjusted  to  be  the  same  as  methane.  The  outer  diameter  of  the  annulus  is  18.0  mm  and 
the  edges  are  thin.  The  annulus  of  combustion  products  from  the  pilot  shields  the  fuel  issuing 
from  the  central  jet  from  contacting  the  surrounding  air  until  a  position  about  4  fuel  jet  diameters 
downstream,  allowing  for  any  unbumed  reactants  in  the  pilot  stream  to  be  consumed.  Sketch  of 
the  piloted  burner  can  be  found  in  Masri  and  Pope’s  paper  (1990).  During  all  the  measurements, 
the  coflow  air  stream  velocity,  ue,  and  the  pilot  burnt  gas  velocity,  uPb,  were  maintained  at  15  and 
24  m/s,  respectively.  The  pilot  burned  gas  velocity,  uPb ,  was  calculated  from  known  reactant  flow 
rates  assuming  a  product  temperature  of  2000  K.  The  flame  characteristics  were  then  changed  by 
varying  the  methane  bulk  jet  velocities  uj. 
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Masri  and  Bilger  (1988)  have  reported  thermocouple  measurements  of  temperatures,  LDA 
measurements  of  velocity  and  turbulence  and  composition  measurements  by  sample  probe. 
Space-  and  time-resolved  measurements  of  temperature,  density,  mixture  fraction,  and  the  mass 
fractions  of  CH4,  O2,  H2O,  H2,  CO,  CO2,  and  N2  have  been  obtained  by  pulsed  spontaneous 
Raman  and  Rayleigh  scattering  in  the  blue  (visibly  soot-free)  regions  of  pure  methane  flames. 
These  data  have  been  reported  in  Masri  et  al.  (1988a,  b,  c).  Table  13  describes  the  various  flames 
investigated  experimentally  and  indicates  the  type  of  measurements  collected  for  each  flame. 
The  overall  equivalence  ratios  for  the  three  fuel  injection  velocities  are  0.043,  0.066,  and  0.088 
from  flame  K,  L,  and  M,  respectively. 


Table  13.  Investigated  Pilot-Stabilized  Diffusion  Flames  of  CH4 


Flame 

Fuel 

Velocity  Uj 

LDA 

Probe 

Raman 

Flame  Characteristics 

K 

27 

V 

V 

All  yellow 

L 

41 

V 

V 

V 

Blue  to  x/D- 30 

M 

55 

V 

V 

V 

Almost  all  blue 

In  this  project,  calculations  are  made  for  methane  flames  at  the  conditions  listed  in  Table  9  (i.e. 
flames  K,  L,  and  M).  The  predictions  are  compared  with  the  relevant  experimental  data,  and  the 
simulation  results  performed  by  Pope  (1990)  using  a  velocity-composition  joint  PDF  transport 
equation  solved  by  the  Monte  Carlo  method.  Here  x  is  the  distance  from  the  burner’s  exit  plane 
and  D  is  the  central  jet  diameter.  In  the  simulations,  the  boundary  conditions  for  all  three 
passages  are  specified  for  the  mean  mixture  fraction,  the  mean  velocities,  and  the  mean  reaction 
progress  variable. 

A  two-step  reaction  mechanism  for  methane  oxidation  was  adopted  for  the  multi-step  model, 

CH4  +  1.5  02  ->  CO  +  2  H20  (4) 

CO  +  0.5  O2  — ^  CO2 

with  Arrhenius  rate  coefficients  used  as, 

q  =  7.943  •  1010  exp(- 15562  IT\CHA\ 573  [02] 008  (5) 

q  =  7.079  •  109  exp(-  5365/T\COf  3[02f 946 

These  rate  parameters  were  reduced  from  the  GRI  3.0  detailed  mechanism  as  discussed  in 
Section  3.4.  Comparisons  with  the  experimental  results  for  velocity,  turbulence,  and  mixing 
fields  for  flames  K,  L,  and  M  are  made  with  the  calculations  from  the  current  multi-step  PDF 
(pdfjms)  model  and  from  the  Monte  Carlo  PDF  (pdf_mc)  model  of  Pope.  Figure  35  shows  the 
simulated  flow  field  by  pdf_ms  for  flame  L.  Fuel  is  injected  from  the  central  hole  at  a  high 
speed,  which  gradually  entrains  the  surrounding  air  to  enter  the  fuel  stream,  causing  the  stream 
velocity  to  decrease  and  the  stream  width  to  expand  down  the  stream.  Fuel  and  air  are  initially 
separated  by  the  burnt  piloted  flame  gas  immediately  surrounding  the  fuel  stream  until  an  axial 
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position  of  5D  is  reached,  where  the  entrained  air  and  fuel  diffuse  toward  each  other  and 
combustion  takes  place  downstream  near  the  fuel  stream  surfaces  as  shown  in  Figure  36  with 
corresponding  rising  temperature. 
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Figure  35.  Calculated  Axial  Velocity  Distribution  for  Flame  L 
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Figure  36.  Calculated  Flame  Temperature  Distribution  for  Flame  L 


Flame  K,  Axial  Velocity  =  27m/s 

Figure  37  shows  plots  of  axial  velocity  (a),  normalized  temperature  (b),  CO  mass  fraction  (c), 
and  O2  mass  fraction  (d)  at  an  axial  location  of  x/d  of  20.  Results  at  an  x/d  of  50  are  shown  in 
Figure  38.  In  general,  agreement  is  very  good.  The  gradient  in  the  02  mass  fraction  is  steeper 
than  measured,  placing  the  flame  location  a  little  farther  radially  out  than  measured.  However, 
CO  predictions  are  relatively  good.  Temperature  predictions  also  match  the  data  well.  Velocity 
matches  are  good,  with  the  center  peak  velocities  being  under  predicted. 
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CO  (’10)  Mass  Fraction  3  Dimensionless  Axial  Velocity 


(a) 


(b) 


(d) 


Figure  37.  Results  for  Diffusion  Flame  K  with  an  Inlet  Axial  Velocity  of27m/s  at  an  x/d  =  20 
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CO  (’10)  Mass  Fraction  Dimensionless  Axial  Velocity 


(a) 


(b) 


(d) 


Figure  38.  Results  for  Diffusion  Flame  K  with  an  Axial  Inlet  Velocity  of27m/s  at  an  x/d  =  50 
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Flame  L,  Axial  Velocity  =  41m/s 

The  results  for  Flame  K,  with  an  axial  velocity  of  41  m/s,  are  plotted  in  Figures  39  and  40. 
Again,  the  overall  agreement  with  the  data  is  relatively  good.  There  is  some  overshoot  in 
temperature,  but  the  results  are  comparable  with  the  more  computationally  expensive  Monte 
Carlo  PDF  method.  The  peak  in  CO  is  over  predicted  and  the  concentration  gradient  on  02  is 
steeper  than  measured.  It  is  likely  that  for  3-step  kinetics  are  insufficient  and  agreement  would 
improve  with  the  application  of  more  complex  mechanisms.  Also,  as  the  jet  velocity  is 
increased,  the  effect  of  local  extinction,  due  to  turbulent  strain,  starts  to  become  important.  This 
effect  is  not  included  in  the  present  model. 


(a) 


(b) 


(c)  (d) 


Figure  39.  Results  for  Diffusion  Flame  L  with  an  Axial  Inlet  Velocity  of  41  m/s  at  an  x/d  =  20 
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Flame  M,  Axial  Velocity  =  55  m/s 

The  results  for  Flame  M,  with  an  axial  velocity  of  55  m/s,  are  plotted  in  Figures  41  and  42. 
Velocity  profile  match  well,  but  the  temperature  and  species  predictions  show  significant 
deviations  from  the  measured  data.  The  same  is  true  for  the  Monte  Carlo  PDF  predictions.  At 
this  jet  velocity,  local  extinction  due  to  turbulent  strain  is  playing  an  important  role.  Since  this 
effect  is  not  accounted  for  in  the  present  model,  the  simulations  show  reactions  taking  place 
earlier  than  measured.  A  turbulent  strain  model  can  be  incorporated  into  the  present  framework 
and  CFDRC  will  look  at  appropriate  options  for  continued  development  of  the  capability  in  the 
future. 


(a) 


(b) 


(c)  (d) 


Figure  41.  Results  for  Diffusion  Flame  M  with  an  Axial  Inlet  Velocity  of  55m/s  at  an  x/d  =  20 
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(a) 


(b) 


(c)  (d) 


Figure  42.  Results  for  Diffusion  Flame  M  with  an  Axial  Inlet  Velocity  of  5  5  m/s  at  an  x/d  =  50 

3.5.2  Premixed  Flame  of  Moreau  et  al. 

The  potential  of  the  multi-step  PDF  model  was  further  studied  by  applying  it  to  a  two- 
dimensional  premixed  combustion  case,  with  the  test  burner  designed  by  Moreau  et  al.  (1976, 
1977).  A  lean-premixed  mixture  of  methane  and  air,  which  has  been  preheated  to  600  K, 
constitutes  the  main  stream  into  the  combustor.  Parallel  to  the  main  stream,  a  pilot  stream  of  hot 
exhaust  gases  at  2000  K  is  injected  for  flame  stabilization.  The  mass  flow  of  the  pilot  stream  is 
approximately  15%  of  the  main  stream.  The  test  section  has  a  length  of  1000  mm.  The  main 
stream  has  a  height  of  80  mm  and  the  pilot  stream  has  a  height  of  20  mm.  The  width  of  the  test 
section  is  100  mm.  The  whole  setup  is  adiabatic,  except  for  some  heat  loss  occurring  at  the  silica 
windows  of  the  test  section. 
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Moreau  et  al.  (1976,  1977)  performed  measurements  for  various  inlet  velocities  and  equivalence 
ratios.  Concentrations  of  chemical  species  were  determined  by  gas  chromatography  of  samples 
taken  with  extractive  probes.  The  most  data  are  available  for  an  equivalence  ratio  of 
approximately  0.82  and  velocities  of  69  m/s  for  the  main  stream  and  121  m/s  for  the  pilot  stream 
close  to  the  inlet  of  the  test  section  (Bohn  and  Lepers,  1999).  Therefore,  these  boundary 
conditions  were  chosen  for  use  in  the  numerical  simulation.  The  experiments  were  performed  at 
atmospheric  pressure.  The  experiments  were  performed  at  atmospheric  pressure,  with  the  flame 
front  observed  to  stabilize  in  the  combustion  chamber  as  shown  in  Figure  43. 


x=0.222  m  x=0.422  m 


Figure  43.  Flame  Front  (Progress  Variance) 


For  this  simulation,  a  two-step  reaction  mechanism  for  methane  oxidation  was  adopted. 

Step  1-CH4+  1.50 2  ->  CO  +  2H20  (29) 

Step  2-CO  +  0.5O2  ->  C02 

with  Arrhenius  rate  coefficients  given  by  Dupont  et  al.  (1993). 

=1.0-1 010  EXP{- 1 20 1 9  IT\CHa  \02  ]  (30) 

q2  =  \.0\0'°EXP(-\2Q>\9IT\CO\02] 

Calculations  using  laminar  chemistry  or  a  PDF  turbulent  combustion  model  with  one-step  global 
mechanism  could  not  predict  a  stable  flame  front  located  inside  the  combustion  chamber  similar 
to  the  experimentally  observed  one.  The  calculated  flame  either  flashes  back  to  the  inlet  or  blows 
out.  The  present  model  predicts  the  flame  front  stabilized  inside  the  combustion  chamber  as 
shown  in  Figure  44. 
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(b)  CO  Mass  Fraction. 
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(c)  CO2  Mass  Fraction 


Figure  44.  Calculated  Distributions  of  Temperature,  CO,  and  CO2  Mass  Fraction 


Comparisons  of  species  and  temperature  profiles  are  shown  in  Figures  45  and  46.  Figure  45 
gives  the  comparison  of  the  calculated  fuel  (CH4)  concentration  distribution  along  the  cross 
section  (y  direction)  with  the  measured  result.  These  fuel  concentration  data  have  been 
normalized  by  the  N2  concentration,  which  is  a  constant  value  in  the  whole  field  as  the  flame  is 
premixed.  At  the  axial  position  of  x=0.222  m,  the  upper  section  is  dominated  with  the  unbumed 
mixture.  The  pilot  flame  ignites  the  premixed  flow,  and  turbulent  combustion  develops  in  the 
mixing  layer  behind  the  board  that  originally  separates  the  unbumed  mixture  with  the  pilot 
flame.  The  current  model  predicts  the  flame  propagation  well  with  good  agreement  of  the  major 
species  concentration  profiles  to  the  measured  data. 
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(a)  x  =  0.422  m 

Figure  45.  Comparison  of  Fuel  Concentration  Profile  Normalized  by  N2  Concentration 


Figure  46  shows  a  comparison  of  calculated  and  measured  flame  temperature  distributions  at 
x=0.422  m.  At  this  location,  the  flame  has  propagated  to  the  upper  section  of  the  main  stream, 
ranging  from  y=0.03  m  to  0.07  m.  Here,  the  present  combustion  model  predicts  an  appropriate 
flame  front  located  around  y=0.05  m.  Some  deviation  between  the  calculated  results  and  the 
measured  data  exist  in  the  high  temperature  region  with  the  model  predicting  a  slightly  faster 
temperature  rise  than  was  measure.  However,  overall  agreement  with  the  experimental  data  is 
good. 
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Figure  46.  Comparison  of  the  Flame  Temperature  Profile 

3.6  Final  Code  Validation 

Several  laboratory  flame  geometries  were  selected  for  final  validation  of  the  newly  developed 
LES  and  multi-step  PDF  combustion  capability.  The  cases  included  a  premixed  flame,  two  swirl 
stabilized  flames,  and  a  jet  diffusion  flame.  The  cases  are  well  established  in  the  literature  and 
are  also  currently  being  used  by  researchers  around  the  world  to  validate  advanced  turbulence- 
combustion  modeling  techniques.  Often  the  flames  are  near  extinction  limits  and  challenging  to 
predict  numerically.  RANS  results  are  typically  poor  and  LES,  with  multi-step  chemistry,  is 
necessary  in  achieving  accurate  results.  Simulations  were  performed  with  the  conventional 
single-step  capability  and  with  the  new  multi-step  PDF  combustion  model.  LES  simulation 
results  with  multi-step  kinetics  were  limited,  due  to  the  inadequacy  of  the  point-solver  approach 
available  in  CFD-ACE+. 

3.6.1  Bluff  Body  Lean-Premixed  Combustor 

Experimental  data  in  a  bluff-body-stabilized  combustor  at  lean  premixed  conditions  were 
obtained  by  Nandula,  Pitz,  and  coworkers  (Pan  et  al.,  1991;  Nandula  et  al.,  1996).  These  laser- 
Doppler  anemometry  (LDA)  measurements  of  velocity  and  a  combination  of  Rayleigh  scattering, 
spontaneous  Raman  scattering,  and  laser-induced  fluorescence  (LIF)  measurements  were 
obtained  to  provide  a  source  of  information  with  which  to  evaluate  combustion  models.  The 
uncertainty  involved  in  the  species  and  temperature  measurements,  based  on  a  single-shot 
uncertainty  analysis,  is  given  in  Table  14. 
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Table  14.  Uncertainty  in  Chemical  Species  Measurements  and  Temperature  Measurements  from 

Nandula  et  al  (1996) 


Scalar  %  Error 


CH4 

3.3 

02 

7.5 

H20 

5.1 

C02 

5.1 

NO 

7.6 

CO 

17.6 

Temperature 

1.0 

The  combustion  chamber,  shown  in  Figure  47,  was  configured  such  that  a  stainless  steel  conical 
bluff  body  of  base  diameter  44.45  mm  and  apex  angle  45°  was  mounted  coaxially  at  the  center  of 
the  combustor  and  served  as  a  flame  holder.  The  burner  enclosure  consisted  of  a  square  cross- 
section  with  filled  round  comers. 


Figure  47.  Schematic  and  Turbulent  Flame  Structure  of  Bluff-Body-Stabilized  Lean  Premixed 

Combustor  (Nandula  et  al.,  1996) 

The  CFD  simulations  use  a  full  3D  geometry  (1015488  cells)  as  shown  in  Figure  48.  The  mass- 
flow  rate  for  the  combustor  was  specified  by  using  the  measured  values  reported  for  the  air  (3960 
slpm)  and  the  CH4  fuel  (244  slpm).  This  corresponded  to  a  CfLj-air  equivalence  ratio  of  0.586. 
The  free  stream  turbulence  at  the  inlet  was  measured  to  be  24%  (based  on  a  15  m/s  mean 
velocity).  Appropriately  scaled  turbulence  levels  were  set  in  the  LES  model  at  the  inlet  for  the 
boundary  condition  of  ksgs.  An  inlet  temperature  of  300  K  and  a  combustor  pressure  of  1  atm 
were  also  specified  in  the  model.  Wall  temperatures  for  the  combustor  and  bluff-body  were  set 
to  500  K. 
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Figure  48.  Structured  Grids  for  3D  Combustor  Geometry 


Table  15  shows  the  multi-step  mechanism  used  for  the  calculations.  The  rates  were  determined 
using  a  fit  to  laminar  ID  flame  calculations  using  the  GRI  3.0  mechanism  as  described  in  Section 
3.2.3. 

Table  15.  Optimized  Arrhenius  Rate  Parameters  Nandula  3-Step  Methane  Oxidation  Mechanism 

(Phi  =  0.586,  T=300  K) 


CH4  +  |o2  =>CO  +  2H20 

/?,  =  4.56  xl09r  01048  exp 

f“: 15273  W]“i”[o2]m“4 

co  +  -o2^co2 

2  z 

R2  =3.76x10 "T0-28  exp 

'  -  8556  VC0]  1.9'^  [0  J  1  0..7 

T  / 

co2^>co+|o2 

R3  =6.023  xl057-a318expf  154261[C02]U34 

\  T  ) 

The  lean  operating  conditions  of  the  Nandula  combustor  make  it  difficult  to  simulate.  The 
RANS  simulations  with  and  without  multi-step  kinetics  did  not  lead  to  a  converged  solution  with 
a  flame.  The  flame  structure  in  the  RANS  simulations  would  gradually  pinch  off  downstream  of 
the  bluff  body  until  the  flame  extinguished.  Figure  49  shows  a  plot  of  the  temperature  profiles  in 
the  RANS  simulations  illustrating  the  pinching  off  and  extinguishing  of  the  flame  downstream  of 
the  bluff  body. 
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Figure  49.  RANS  Simulation  Results  Prior  to  Blowout  Showing  the  Extinguishing  of  the  Flame 

Downstream  of  the  Bluff  Body 

By  capturing  the  turbulent  flow  structures  that  mix  hot  products  and  reactants  together,  LES 
simulations  are  capable  of  predicting  a  stable  flame.  At  the  very  lean  operating  conditions  of  this 
validation  case,  LES  is  essential  to  capture  realistic  flame  behavior.  The  simulation  is  also  very 
sensitive  to  the  kinetics  used.  Due  to  bug  fixes  and  software  problems,  the  simulation  was 
restarted  multiple  times.  Figures  50  through  52  shows  snapshots  from  the  LES  calculations  for 
temperature,  axial  velocity,  and  CO  mass  fraction.  Averaged  quantities  are  not  yet  available  for 
the  simulation. 


Figure  50.  Snapshot  of  Temperature  in  the  Nandula  Bluff-Body  Combustor  LES  Simulation 

Using  Multi-Step  Chemistry 
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Figure  51.  Snapshot  of  Axial  Velocity  in  the  Nandula  Bluff-Body  Combustor  LES  Simulation 

Using  Multi-Step  Chemistry 
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Figure  52.  Snapshot  of  CO  Mass  Fraction  in  the  Nandula  Bluff-Body  Combustor  LES 

Simulation  Using  Multi-Step  Chemistry 


Comparison  of  the  predicted  temperature,  CO,  CO2,  H2O,  and  O2  are  shown  in  the  following 
plots  (Figures  53  through  56)  at  axial  locations  with  an  x/d  of  0.1,  0.6,  and  1.0.  Only 
instantaneous  LES  results  are  available  at  this  time.  As  the  results  are  averaged,  the  peaks  in  the 
profiles  will  be  reduced  as  the  fluctuations  in  the  flame  surface  are  accounted  for.  Overall  the 
radial  flame  location  is  predicted  well.  As  expected,  the  gradients  in  the  instantaneous  LES  plots 
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are  much  steeper  than  those  observed  in  the  averaged  data.  The  fluctuations  captures  within  the 
LES  simulation  will  tend  to  smooth  out  gradients  and  decrease  the  peak  magnitudes  as  the  results 
are  averaged.  CO  concentrations  downstream  of  the  bluff  body  are  significantly  lower  (about 
80%)  than  the  measured  data. 


Figure  53.  Plots  of  Instantaneous  LES  Results  Compared  to  the  Averaged  Measured  Data  for 
CO,  T,  CO2,  and  H2O  Downstream  of  the  Bluff  Body  at  an  x/d  of  0.1 
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CO  (Mass  Fraction)  T  (K) 


Figure  54.  Plots  of  Instantaneous  LES  Results  Compared  to  the  Averaged  Measured  Data  for 
CO,  T,  CO  2,  and  H2O  Downstream  of  the  Bluff  Body  at  an  x/d  of  0.6 


68 


8503/8 


c 

o 

as 

o 

co 


</> 

w 

(0 

2 


r/d 


r/d 


Figure  55.  Plots  of  Instantaneous  LES  Results  Compared  to  the  Averaged  Measured  Data  for 
CO,  T,  CO  2,  and  HiO  Downstream  of  the  Bluff  Body  at  an  x/d  of  1.0 


3.6.2  SMA2  -  Swirl  Stabilized  Methane  Air  Flame 

Swirl  stabilized  diffusion  flames  are  often  observed  in  practical  combustors.  The  new  CFD 
capabilities  developed  in  this  program  must  be  able  to  resolve  the  characteristics  of  these  flames 
correctly.  As  part  of  the  TNF  workshop,  detailed  experiments  have  been  performed  for  swirl 
stabilized  diffusion  flames.  A  comprehensive  set  of  validation  data  has  been  collected  that  can  be 
used  to  validate  turbulence  and  chemistry  models.  For  this  reason,  the  SMA2  flame  was  selected 
as  a  validation  case.  The  SMLA2  flame  is  the  high  swirl  diffusion  flame  extensively  documented 
in  the  TNF  workshop.  The  following  sections  briefly  review  the  burner  details  presented  in  the 
last  quarterly  report.  The  CFD  setup  and  results  from  the  RANS  calculations  are  presented  in  the 
next  section. 

Experimental  Setup:  The  swirl  burner  for  the  SMA2  case  is  centrally  located  in  an  air  co-flow. 
The  co-flow  is  confined  in  a  square  duct  with  sides  130  mm  in  length.  Air  flows  through  the 
swirler  annulus,  which  is  supplied  through  the  main  inlet  and/or  the  tangential  injectors  as  shown 
in  Figure  56.  Figure  57  shows  a  picture  of  the  SMA2  flame  [TNF  workshop,  2002].  The  fuel  (2 
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parts  air  and  1  part  methane  by  volume)  issues  from  the  central  jet,  located  in  the  center  of  the 
ceramic  bluff-body  face.  The  jet  is  fully  developed  at  the  burner  exit.  The  SMA2  flame 
characteristics  are  summarized  in  Table  16. 


Figure  56.  SMA2  Burner  Schematic 
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Figure  57.  The  SMA2  Flame 


Table  16.  SMA2  Operating  Conditions. 


Fuel 

mixture 

Axial 

annulus 

velocity 

(m/s) 

Tangential 

annulus 

velocity 

(m/s) 

Fuel  jet 
velocity 
(m/s) 

Reynolds 
number  at 
annulus  exit 

Reynolds 
number 
of  fuel  jet 

Swirl 

number 

Methane 

-Air 

16.3 

25.9 

66.3 

32400 

15400 

1.59 

The  experimental  data  set  for  this  case  consists  of  radial  profiles  of  axial  and  tangential  velocity, 
temperature  and  mass  fraction  of  OH,  CO,  CO2,  H2O  at  many  locations  downstream  of  the 
burner  exit.  This  set  of  data  results  in  a  comprehensive  test  of  the  CFD  capabilities  to  predict  the 
characteristics  of  a  swirl-stabilized  diffusion  flame. 

CFD-Model  Description:  The  experimental  data  for  the  SMA2  flame  was  collected  starting 
from  6.8  mm  downstream  of  the  burner  exit.  To  get  the  boundary  conditions  for  the  3D  case  at 
the  burner  exit  for  the  jet,  annulus  and  co-flow,  a  2D  axi-symmetric  non-reacting  RANS 
simulation  was  performed.  The  2D  axi-symmetric  grid  is  shown  in  Figure  58. 
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Figure  58.  2D  Axi-Symmetric  Grid  for  Non-Reacting  SMA2  Case 


The  swirler  and  fuel  inlets  were  made  sufficiently  long  to  allow  the  flow  to  fully  develop  at  the 
burner  exit.  The  boundary  conditions  at  the  swirler  and  jet  inlets  were  varied  till  the  velocity 
profile  from  the  CFD  simulations  matched  the  velocity  profile  measured  experimentally  at  6.8 
mm  downstream  of  the  burner  exit.  This  provided  a  means  to  estimate  the  boundary  conditions 
for  the  swirler  and  jet  inlets  in  the  3D  reacting  RANS  case.  The  combustor  walls  were  given 
adiabatic  boundary  conditions  and  a  constant  pressure  outlet  was  used  for  the  exit  boundary.  The 
standard  k-epsilon  turbulence  model  was  used  for  turbulence  closure.  The  comparison  between 
the  experimental  measurements  and  2D  axi-symmetric  CFD  RANS  simulation  for  the  non¬ 
reacting  case  is  shown  in  Figure  59. 


Figure  59.  Axial  and  Tangential  Velocity  Comparisons  Between  Experimental  and  2D  Axi- 
Symmetric  CFD  Results  at  6.8  mm  Downstream  of  Burner  Exit 


The  results  shown  in  Figure  59  indicated  that  the  right  boundary  conditions  for  the  burner  exit 
were  now  available  to  do  a  full  3D  simulation.  Excellent  agreement  in  axial  velocity  was 
achieved.  Agreement  to  the  swirl  velocity  was  similar  to  what  has  been  report  by  other 
researchers  (Masri,  2002).  A  3D  steady  state  RANS  non-reacting  case  was  then  set  up  to  get  the 
initial  conditions  for  the  reacting  case.  The  3D  grid  is  shown  in  Figure  60.  The  computational 
grid  had  a  total  of  1,263,600  cells. 
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Figure  60.  3D  RANS  Grid,  (a)  Top  Left:  Front  View,  (b)  Top  Right:  Close-Up  of  Front  View  and 

(c)  Bottom:  Side  View 


The  standard  k-epsilon  model  was  used  for  turbulence  closure  with  Schmidt  and  Prandtl  numbers 
of  0.7.  The  inlet  boundary  conditions  for  the  swirler  annulus  and  the  jet  inlet  were  taken  from  the 
2D  axi-symmetric  case  previously  discussed.  Central  differencing  was  used  for  all  spatial 
discretizations.  Figures  61  shows  the  axial  velocity  contours  of  the  RANS  result. 


8  ?  8  S 

Figure  61.  RANS  Cold-Flow  Axial  Velocity 


Figures  62-64  show  the  comparisons  between  the  experimental  measurements  and  CFD 
predictions  of  air  velocity  at  different  axial  locations  downstream  of  the  burner  exit. 
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Figure  62.  Axial  and  Tangential  Velocity  Comparisons  Between  Experimental  and  3D  CFD 
RANS  Results  at  6.8  mm  Downstream  of  Burner  Exit 


Figure  63.  Axial  and  Tangential  Velocity  Comparisons  Between  Experimental  and  3D  CFD 
RANS  Results  at  30.0  mm  Downstream  of  Burner  Exit 


y(mm)  y(mm) 


Figure  64.  Axial  and  Tangential  Velocity  Comparisons  Between  Experimental  and  3D  CFD 
RANS  Results  at  100.0  mm  Downstream  of  Burner  Exit 

The  CFD  results  show  good  comparison  to  data  for  the  non-reacting  case.  Some  detailed  features 
adjacent  to  the  inlets  are  not  well  resolved,  due  to  the  coarser  grid  spacing  of  the  3D  mesh.  With 
the  non-reacting  case  converged  and  matching  well  with  experimental  data,  a  reacting  flow  case 
was  performed. 
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The  reacting  case  used  the  same  grid  as  was  used  by  the  non-reacting  case  above.  To  start  with, 
the  k-epsilon  model  was  used  for  the  turbulence  closure  and  the  single  step  assumed  PDF  model 
was  used  for  modeling  the  effects  of  turbulence-chemistry  interaction.  The  3D  reacting  RANS 
case  was  initialized  by  the  flow  field  obtained  from  the  non-reacting  case.  Central  differencing 
was  used  for  all  spatial  discretizations.  Contour  plots  of  axial  velocity  and  temperature  are  shown 
in  Figures  65  and  66.  Comparison  to  data  at  various  axial  locations  of  the  3D  steady  state, 
reacting-flow,  and  RANS  simulations  are  shown  in  Figure  67. 


Figure  65.  RANS  Axial  Velocity  -  Reacting  Flow 


Figure  66.  RANS  Temperature  -  Reacting  Flow 
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Figure  67.  Axial  and  Tangential  Velocity  Comparisons  (Top)  and  Temperature  Comparisons 
( Bottom )  Between  Experimental  and  3D  CFD  Results  at  10.0  mm  Downstream  of  Burner  Exit 

Figure  68  shows  the  comparison  of  the  radial  profiles  of  velocity  components  and  temperature  of 
the  flow-field  between  the  CFD  results  and  the  experimentally  collected  data  at  10  mm 
downstream  of  the  burner  exit.  The  axial  velocity  predictions  compare  well  with  experimental 
measurements.  However,  the  RANS  simulation  is  not  able  to  capture  the  two  peaks  in  the 
tangential  velocity  as  noted  in  the  experiments.  Traditionally  RANS  simulations  have  not 
performed  very  well  for  highly  swirling  and  reacting  flows.  The  temperature  profile  shape  and 
magnitudes  are  also  not  well  captured  by  the  RANS  simulation.  The  peak  temperature  is  under¬ 
predicted  by  350K,  while  the  hot  recirculation  zone  near  the  centerline  is  not  captured  quite  as 
well.  These  results  are  similar  to  RANS  results  for  this  flame  reported  by  Masri  et  al.  (2002). 

Figure  69  shows  the  comparisons  of  the  CFD  predictions  with  the  experimental  measurements 
30  mm  downstream  of  the  burner  exit.  The  axial  and  tangential  velocities  shapes  are  also  not 
well  captured  by  the  RANS  predictions.  The  velocity  profiles  predicted  by  the  CFD  calculations 
have  been  substantially  dissipated  near  sharp  changes  in  the  experimental  observed  profiles. 
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Figure  68.  Axial  and  Tangential  Velocity  Comparisons  (Top)  and  Temperature  Comparisons 
( Bottom )  Between  Experimental  and  3D  CFD  Results  at  30.0  mm  Downstream  of  Burner  Exit 


y  (mm) 

Figure  69.  Axial  and  Tangential  Velocity  Comparisons  (Top)  and  Temperature  Comparisons 
(Bottom)  Between  Experimental  and  3D  CFD  Results  at  100.0  mm  Downstream  of  Burner  Exit 
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Figure  70  also  shows  the  decay  in  the  CFD  predictions  of  the  RANS  simulation  compared  to  the 
experimentally  observed  profiles.  Such  a  behavior  is  representative  of  RANS  calculations  for 
highly  swirling  and  reacting  flows.  It  is  expected  that  the  LES  results  will  be  able  to  capture  the 
relevant  flow  features  much  better  than  RANS. 

A  3D  LES  calculation  using  the  assumed  PDF  model  was  performed  for  the  SMA2  case.  The 
LDKM  model  was  used  for  the  subgrid  closure.  The  grid  used  for  the  LES  case  is  the  same  as 
RANS  cases  above.  The  Crank-Nicolson  scheme  with  a  blending  factor  of  0.6  is  used  for  time 
integration  while  the  central  differencing  scheme  is  employed  for  spatial  discretizations.  Figure 
22  shows  a  snapshot  of  the  LES  temperature  contours  for  the  SMA2  case. 


Figure  70.  Snapshot  of  Temperature  for  the  LES  SMA2  Case 


Comparison  to  data  for  two  downstream  axial  locations  of  the  3D  reacting  LES  simulation  is 
shown  in  Figure  71-72. 
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Figure  7L  Axial  Velocity  ( Top  Left),  Tangential  Velocity  (Top  Right),  CO2  Mass  Fraction 
( Bottom  Left )  and  Temperature  ( Bottom  Right)  Comparisons  Between  Experimental  and  3D 
CFD  Results  at  10.0  mm  Downstream  of  Burner  Exit 
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Figure  72.  Axial  Velocity  ( Top  Left),  Tangential  Velocity  (Top  Right),  CO 2  Mass  Fraction 
(Bottom  Left)  and  Temperature  (Bottom  Right)  Comparisons  Between  Experimental  and  3D 
CFD  Results  at  10.0  mm  Downstream  of  Burner  Exit 


The  LES  results  shows  signs  of  improvement  over  the  RANS  results.  The  axial  velocity  profiles 
are  almost  perfectly  captured,  while  the  tangential  velocity  profiles  seem  to  miss  the  first  peak  as 
observed  in  the  experiments.  The  species  predictions,  although  better  than  RANS,  are  not 
captured  correctly.  The  temperature  profile  with  LES  is  flatter  and  more  consistent  with  the 
observed  profile. 

2D  simulations  were  performed  using  both  the  3  and  5  step  mechanisms  shown  in  Section  3.2. 
Also,  results  using  the  conventional  single-step  chemistry  are  shown  for  comparison  purposes. 
Axial  and  tangential  velocity  results  are  shown  in  Figures  73  and  74,  plotted  at  axial  locations  of 
6.8,  30,  and  100  mm  downstream  of  the  swirler  exit.  Initial  results  close  to  the  boundary  are 
good,  but  worsen  farther  downstream.  Typically  RANS  has  a  difficult  time  with  swirling  flows, 
and  these  difficulties  are  highlighted  in  the  results.  Swirling  flow  velocities  near  the  centerline 
underpredicted.  This  has  also  been  shown  by  other  researchers  for  RANS  and  LES  is  expected 
to  result  in  improved  velocity  predictions.  As  shown  in  Figure  73c,  chemistry  also  has  a  strong 
influence  on  the  centerline  velocity.  Overall  results  with  improved  chemistry  result  in  improved 
velocity  predictions. 
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Figure  73.  SMA2  Flame  Axial  (u-)  Velocity  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x  =  6.8mm  (b)  x  =  30mm  (c)  =  100mm 
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(a) 


(b) 


(c) 


Figure  74.  SMA2  Flame  Tangential  (w-)  Velocity  Comparisons  at  Various  Axial  (x-) 
Locations  Downstream  of  the  Burner  (a)  x  =  6.8mm  (b)  x  =  30mm  (c)  =  100mm 


Temperature,  CO2,  and  CO  results  are  plotted  in  Figures  75  through  77.  Near  the  swirler,  only 
minor  differences  with  different  chemical  mechanisms  are  noted,  with  temperatures  reaching  the 
experimental  values  only  near  the  centerline.  Downstream,  peak  temperatures  are  much  more 
realistic  with  both  the  3  and  5  step  mechanisms.  Initial  peak  CO  values  are  higher  than 
measured.  However,  the  multi-step  mechanisms  do  a  much  better  job  of  predicting  CO 
concentrations  at  the  100  mm  location.  Here,  the  single-step  mechanism  results  in  essentially 
zero  CO,  while  the  multi-step  mechanism  results  show  high  CO  values  on  the  centerline.  The 
results  show  that  multi-step  mechanisms  provide  improved  predictions,  but  there  are  still 
significant  discrepancies  between  the  predictions  and  the  measured  results.  It  is  expected  that 
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Temperature  (K)  w  Temperature  (K) 


the  application  of  LES  will  result  in  improved  mixing  and  velocity  predictions,  and  therefore  will 
provide  improved  temperature  and  concentration  results. 
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Figure  75.  SMA2  Flame  Temperature  Comparisons  at  Various  Axial  ( x -)  Locations 
Downstream  of  the  Burner,  (a)  x  —  10mm  (b)  x  =  30mm  (c)  =  100mm 
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Figure  76.  SMA2  Flame  CO  Mass  Fraction  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x  =  10mm  (b)  x  =  30mm  (c)  =  100mm 
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Figure  77.  SMA2  Flame  CO2  Mass  Fraction  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x  =  10mm  (b)  x  =  30mm  (c)  =  100mm 


The  3D  model  of  the  SMA2  geometry  with  multi-step  chemistry  was  attempted  but  a  converged 
solution  was  not  obtained.  The  unswirled  center  jet  make  obtaining  converged  3D  simulations 
difficult  and  the  difficulty  is  significantly  increased  using  the  multi-step  PDF  model.  As 
previously  mentioned,  the  solvers  available  in  CFD-ACE+  have  proven  to  be  inadequate  for 
solving  the  coupled  species  and  energy  equations.  The  slow  propagation  of  species  information 
relative  to  the  enthalpy  results  in  unrealistically  low  temperatures  in  the  flow  domain,  preventing 
convergence. 
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3.6.3  Sandia-D  -  Piloted  let  Flame 

The  third  validation  case  selected  is  the  Sandia-D  flame.  This  flame  has  also  been  well 
characterized  as  part  of  the  TNF  workshop  and  is  widely  used  as  a  validation  case  for  diffusion 
flames.  As  opposed  to  the  swirl-stabilized  SMA2  flame  discussed  above,  the  Sandia-D  flame  is 
a  pilot-stabilized  flame.  Figure  78  shows  a  close-up  picture  of  the  Sandia-D  flame. 

The  Sandia-D  flame  provides  a  good  validation  case  to  test  the  turbulence-chemistry  interaction 
model  in  conditions  often  observed  in  practical  combustors.  A  brief  summary  of  the 
experimental  setup  is  presented  below. 


Figure  78.  The  Sandia-D  Flame 


The  piloted,  axi-symmetric  burner  has  a  main  jet  diameter  of  7.2  mm  and  a  pipe  length 
exceeding  40  times  the  diameter,  thus  ensuring  a  fully  developed  flow  at  the  pipe  exit.  The  jet 
composition  is  25%  methane  and  75%  air  by  volume.  The  jet  exit  Reynolds  number  is  22,400. 
The  diameter  of  the  annular  pilot  is  18.2  mm  and  consists  of  72  tiny  premixed  jets  creating  a 
homogenous  distribution  of  temperature  and  gas  composition.  A  mixture  of  acetylene,  hydrogen, 
air,  carbon  dioxide  and  nitrogen  is  chosen  with  the  same  enthalpy  and  equilibrium  composition 
as  methane/air  at  0.77  equivalence  ratio  for  the  pilot.  The  burner  is  placed  concentrically  to  the 
exit  nozzle  of  a  wind  tunnel  that  provides  a  laminar  annular  co-flow  of  filtered  air  at  a  speed  of 
approximately  1.0  m/s.  The  jet  exit  is  placed  20  mm  above  the  exit  of  the  co-flow  nozzle.  All  gas 
flows  before  combustion  are  at  room  temperature. 

CFD-Model  Description:  The  experimental  data  for  the  Sandia-D  case  was  collected  starting 
from  1mm  downstream  of  the  burner  exit,  thus  providing  a  good  boundary  condition  to  start  the 
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calculations.  Experimental  profiles  of  the  axial  and  radial  velocities  were  specified  at  the  fuel 
and  air  inlets.  A  constant  axial  velocity  of  0.9  m/s  was  specified  for  the  co-flow.  The  3D  grid  for 
the  CFD  calculation  is  shown  in  Figure  79. 


Figure  79.  3D  RANS  Grid,  (a)  Top  Left:  Front  View,  (b)  Top  Right:  Close-Up  of  Front  View  and 

(c)  Bottom:  Side  View 


The  RANS  simulation  used  the  k-e  turbulence  model  and  the  single  step  assumed  PDF  model  for 
turbulence-chemistry  closure.  Upwind  differencing  was  used  for  all  spatial  discretizations. 
Contour  plots  of  axial  velocity  and  temperature  are  shown  in  Figures  80  and  81.  Comparison  to 
data  at  various  axial  locations  of  the  3D  steady  state,  reacting-flow,  RANS  simulation  are  shown 
in  Figures  82-86. 
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Figure  82.  Axial  Velocity  (Top  Left),  Temperature  (Top  Right),  CH4  Mass  Fraction  (Bottom 
Left)  and  CO  Mass  Fraction  (Bottom  Right)  Comparisons  Between  Experimental  and  3D  CFD 

Results  at  xJd  =  1  Downstream  of  Burner  Exit 

As  shown  in  Figure  82,  the  RANS  results  show  excellent  predictions  for  velocity,  temperature 
and  species  just  downstream  of  the  burner  exit.  The  CO  mass  fraction  is  over  predicted,  which 
will  be  corrected  when  the  multi-step  chemistry  closure  is  used. 

As  we  move  downstream  (Figure  83  and  84),  the  velocity  predictions  still  show  a  great  match 
with  data,  but  the  temperature  predictions  deteriorate.  The  reason  for  this  discrepancy  can  be 
attributed  to  the  high  resolution  required  where  the  pilot  holds  the  flame.  RANS  simulations 
cannot  provide  the  rapid  mixing  and  burning  predictions  in  this  zone.  LES  provides  the 
necessary  capabilities  to  resolve  this  region.  Due  to  the  temperature  under-predictions,  the 
species  are  also  not  well  predicted  in  this  region. 
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Figure  83.  Axial  Velocity  (Top  Left),  Temperature  ( Top  Right),  CH4  Mass  Fraction  (Bottom 
Left)  and  CO  Mass  Fraction  (Bottom  Right)  Comparisons  Between  Experimental  and  3D  CFD 

Results  at  x/d  =  3  Downstream  of  Burner  Exit 
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Figure  84.  Axial  Velocity  ( Top  Left),  Temperature  ( Top  Right),  CH4  Mass  Fraction  (Bottom 
Left )  and  CO  Mass  Fraction  ( Bottom  Right)  Comparisons  Between  Experimental  and  3D  CFD 

Results  at  x/d  =  15  Downstream  of  Burner  Exit 


Further  downstream,  RANS  predictions  again  improve  and  velocity,  temperature  and  species  are 
well  predicted  (Figures  85  and  86).  Downstream  of  this  location,  methane  is  all  consumed  and 
thus  only  CO  predictions  are  shown. 
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Figure  85.  Axial  Velocity  (Top  Left),  Temperature  ( Top  Right),  and  CO  Mass  Fraction  ( Bottom ) 
Comparisons  Between  Experimental  and  3D  CFD  Results  at  x/d  =  30  Downstream  of  Burner 

Exit 
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Figure  86.  Axial  Velocity  ( Top  Left),  Temperature  ( Top  Right),  and  CO  Mass  Fraction  ( Bottom ) 
Comparisons  Between  Experimental  and  3D  CFD  Results  at  x/d  =  45  Downstream  of  Burner 

Exit 


A  number  of  RNG  k-e  turbulence  model  based  RANS  simulations,  both  2D  axi-symmetric  and 
full  3D,  were  performed  on  the  Sandia-D  flame  using  multi-step  chemistry.  The  simulations 
were  carried  out  using  both  the  1-step  mechanism  as  well  as  the  multi-step  assumed  PDF  model 
for  turbulence-chemistry  closure.  Both  the  3-step  and  5-step  mechanisms  described  in  Section 
3.4  were  used.  Upwind  differencing  was  used  for  all  spatial  discretizations.  The  2D  axi- 
symmetric  and  the  3D  RANS  data  were  compared  with  the  experimental  data  at  various  locations 
along  the  axis  of  the  flame.  The  grid  was  also  modified  and  improved  from  what  was  used  in  the 
simulation  results  shown  previously.  The  grid  features  fine  grid  resolution  near  the  centerline  to 
capture  the  small  center  fuel  jet  and  is  adequate  for  both  RANS  and  LES  simulations.  The  grid 
is  shown  in  Figure  87. 
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Figure  87.  Grid  for  Both  RANS  and  LES  Simulations  of  the  Sandia-D  Flame 


Plots  of  temperature  and  axial  velocity  for  the  multi-step  RANS  simulation  are  shown  in  Figures 
88  and  89.  Plots  of  CO  and  CO2  mass  fractions  are  shown  in  Figure  90  and  91.  Since  the 
measurements  were  taken  in  unconfined  flow,  the  outer  boundaries  of  the  model  were  set  to 
constant  pressure  exit  boundaries.  Table  17  shows  the  5-step  mechamism  used  for  the 
simulation.  The  mechanism  was  generated  from  the  GRI  3.0  detailed  mechanism  using  a  series 
of  PSR/PFR  calculations  over  an  equivalence  ratio  range  of  0.6  to  1.4  as  described  in  Section  3.3 


Table  17.  Five-Step  Methane  Mechanism  Accounting  for  CO  and  H2 


Reaction 

ln(A) 

Ea/R 

TK1 

P 

Reactant  Order 

2  CH4  +  02  =>  2  CO  +  4  H2 

34.53 

35872 

2.618  CH4 

1.109  02 

1.572 

2  H2  +  02  =>  2  H20 

34.53 

8193 

-1.421  02 

1.363  H2 

1.063 

2  CO  +  02  =>  2  C02 

34.53 

13574 

-0.098  02 

1.601  CO 

1.318 

2  H20  =>  2  H2  +  02 

12.01 

50000 

0.093  H20 

1.242 

2  C02  =>  2  CO  +  02 

34.53 

48202 

-0.177  C02 

1.507 
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Figure  88.  RANS  Temperature  Profiles  for  the  Sandia-D  Flame  Using  Multi-Step  Chemistry 


Figure  89.  RANS  Axial  Velocity  Profiles  for  the  Sandia-D  Flame  Using  Multi-Step  Chemistry 
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Figure  90.  RANS  CO  Mass  Fraction  Profiles  for  the  Sandia-D  Flame  Using  Multi-Step 

Chemistry 


Figure  91.  RANS  CO 2  Mass  Fraction  Profiles  for  the  Sandia-D  Flame  Using 

Chemistry 


C02 


Multi-Step 


Figure  92  shows  a  comparison  between  the  numerical  and  the  experimental  axial  velocity  data  at 
four  axial  locations.  It  can  be  seen  that  among  the  2D  axisymmetric  cases,  the  5-step  mechanism 
case  compares  best  with  experiments.  As  is  to  be  expected,  the  3D  5-step  mechanism  case  comes 
closest  to  the  experimental  data  among  all  the  numerical  cases.  In  Figure  92(a),  the  3D  5-step 
mechanism  case  is  able  to  resolve,  to  some  extent,  the  sharp  gradients  in  the  axial  velocity 
around  r/D  =  1.  However,  in  Figure  92(c),  at  x/D  =  30,  the  3D  RANS  case  over  predicts  the  peak 
axial  velocity  at  r/D  =  0,  even  though  it  compares  very  well  with  the  experimental  data  at  larger 
radial  distances. 
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(a)  (b) 


Figure  92.  Sandia-D  Flame  Axial  (u-)  Velocity  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x/d  =  1  (b)  x/d  =  15  (c)  x/d  =  30  (d)  x/d  =  45 


A  comparison  of  flame  temperatures  at  various  axial  locations  is  shown  in  Figure  93.  Again,  the 
3D  5-step  RANS  case  does  the  best.  The  improvement  in  predictions  compared  to  the 
corresponding  2D  axisymmetric  5-step  case  is  quite  significant.  At  all  the  four  axial  locations, 
the  3D  5-step  case  shows  remarkably  good  agreement  with  the  experiments,  especially  the  peak 
temperatures.  Even  at  far  downstream  locations  such  as  at  x/D  =  45,  shown  in  Figure  93(d),  the 
3D  case  does  a  much  better  job  than  all  the  2D  cases.  As  is  to  be  expected,  2D  1-step  RANS  case 
performs  rather  poorly  in  predicting  both  the  magnitude  and  the  radial  locations  of  the  peak 
temperatures  at  all  the  axial  locations  shown. 
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Figure  93.  Sandia-D  Flame  Temperature  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner,  (a)  x/d  -  1  ( b )  x/d  =  15  (c)  x/d  =  30  (d)  x/d  =  45 


The  CH4  species  mass  fractions  are  shown  in  Figure  94.  At  the  first  two  axial  locations,  x/D  =  1 
and  x/D  =15,  again  the  3D  case  compares  well  with  the  experiments,  especially  in  predicting  the 
peak  CH4  mass  fractions.  However,  at  subsequent  downstream  locations,  the  comparison 
between  the  3D  and  the  experimental  data  deteriorates.  Both  the  2D  and  the  3D  5-step  chemistry 
cases  indicate  a  near-complete  burnout  of  CH4  at  the  two  farther  downstream  axial  locations.  The 
CO  mass  fraction  profiles  are  shown  in  Figure  95.  Here,  none  of  the  2D  or  the  3D  cases  shows 
good  comparison  with  experiments.  However,  in  general,  the  multi-step  chemistry  cases  better 
capture  the  trends  in  the  CO  concentration  profiles. 
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Finally,  the  numerical  and  experimental  data  for  the  CO2  mass  fractions  are  compared  in  Figure 
96.  The  multi-step  chemistry  cases  show  an  improvement  in  predictions.  In  particular,  the  2D 
and  the  3D  5-step  mechanism  cases  accurately  capture  the  radial  locations  of  the  peak  CO2 
concentrations,  even  though  the  peak  magnitudes  are  not  captured  well.  The  1-step  mechanism 
case,  as  usual,  performs  the  worst. 
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Figure  94.  Sandia-D  Flame  CH4  Mass  Fraction  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x/d  =  1  (b)  x/d  =  15  (c)  x/d  =  30 
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Figure  95.  Sandia-D  Flame  CO  Mass  Fraction  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x/d  =  1  (b)  xJd  =  15  (c)  x/d  =  30  (d)  x/d  =  45 
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Figure  96.  Sandia-D  Flame  CO 2  Mass  Fraction  Comparisons  at  Various  Axial  (x-)  Locations 
Downstream  of  the  Burner  (a)  x/d  =  1  (b)  x/d  =  15  (c)  x/d  =  30  (d)  x/d  =  45 


LES  simulations  of  the  Sandia-D  flame  are  currently  ongoing,  but  have  not  proceeded  far  enough 
for  a  meaningful  comparison  to  the  data.  Preliminary  results  from  this  simulation  are  shown  in 
Figures  97  through  99.  The  flame  is  still  significantly  longer  than  what  was  observed  in  the 
RANS  simulations,  but  as  the  simulations  proceeds  toward  a  limit-cycle,  the  increase  in  the 
turbulence  will  tend  to  shorten  and  spread  out  the  flame.  Significant  vortices  are  seen 
developing  in  the  jet,  as  illustrated  by  the  fluctuations  in  the  radial  velocity  components. 
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Figure  97.  Preliminary  Snapshot  of  Temperature  from  the  Sandia-D  LES  Combustor 

Simulations 
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Figure  99.  Preliminary  Snapshot  of  Radial  Velocity  from  the  Sandia-D  LES  Combustor 

Simulations 


3.6.4  TECFLAM  Combustor  -  Swirling  Diffusion  Flame 

The  fourth  validation  case  selected  is  the  TECFLAM  combustor.  This  case  was  selected  to 
provide  a  more  in-depth  evaluation  of  the  chemistry  model  for  NOx  predictions.  The  TECFLAM 
consortium  in  Germany  has  collected  a  comprehensive  set  of  experimental  data  for  a  swirl- 
stabilized  combustor  for  validation  of  CFD  turbulence  and  chemistry  models.  The  TECFLAM 
combustor  provides  detailed  measurements  on  NOx  formation  in  swirl-stabilized  combustors  and 
thus  forms  a  good  validation  case  for  turbulence-chemistry  interaction  models  for  NOx 
predictions.  A  schematic  of  the  TECFLAM  combustor  is  shown  in  Figure  100. 
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Figure  100.  Schematic  of  the  TECFLAM  Combustor  (Meier  et  al.,  2000) 
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Natural  gas  and  air  are  supplied  through  the  annular  nozzles.  The  inside  diameter  for  fuel  is  20 
mm  while  the  outside  diameter  is  26  mm.  The  inside  diameter  for  air  is  30  mm  while  the  outside 
diameter  is  60  mm.  The  fuel  contains  96.4-98.2%  methane  by  volume  and  minor  admixtures  of 
nitrogen,  carbon  dioxide  and  higher  hydrocarbons.  The  flame  has  a  thermal  power  of  150  kW, 
equivalence  ratio  of  0.833  and  a  swirl  number  of  0.9.  The  Reynolds  number  at  the  air  inlet  is 
42,900  while  at  fuel  inlet  is  7900.  The  water  cooled-combustion  chamber  (Twater_iii  =  65°C, 
TWater_out=  75°C)  has  an  inner  diameter  of  500  mm  and  an  axial  length  of  1200  mm.  Figure  101 
shows  a  picture  of  the  TECFLAM  combustor  and  flame. 


CFD-Model  Description:  The  experimental  data  for  the  TECFLAM  case  was  collected  starting 
from  1mm  downstream  of  the  burner  exit,  thus  providing  a  good  boundary  condition  to  start  the 
calculations.  Experimental  profiles  of  the  axial  and  tangential  velocities  were  specified  at  the  air 
inlets,  while  a  bulk  velocity  of  21  m/s  was  specified  for  the  fuel  inlet.  A  constant  heat  flux  of  50 
kW  was  specified  for  the  liners.  The  3D  grid  for  the  CFD  calculation  is  shown  in  Figure  102. 
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Figure  102.  3D  RANS  Grid,  (a)  Top  Left:  Angled  View,  (b)  Top  Right:  Side  View  and  (c) 

Bottom:  Front  View 


CFD  Results:  The  initial  RANS  predictions  with  single-step  chemistry  show  decent  velocity 
predictions,  but  the  temperature  predictions  show  a  poor  match.  The  lack  of  appropriate  heat 
flux  boundary  conditions  on  the  combustor  liners  is  probably  responsible  for  this  discrepancy. 
Figures  103-107  show  the  results  from  the  3D  reacting  steady  state  RANS  simulation.  Although 
the  axial  velocity  shows  decent  comparison  with  data  at  all  axial  locations  compared,  the 
tangential  velocity  and  temperatures  are  very  poorly  predicted.  RANS  simulations  are  known  to 
perform  poorly  in  predicting  the  swirl  velocity  component.  The  temperature  profiles  are  also 
mismatched  due  to  the  incorrect  flow  pattern  that  sets  up  because  of  the  flow  field  prediction 
errors  in  RANS. 
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Figure  103.  Comparison  of  3D  RANS  with  Experimental  Data  at  10  mm  Downstream  of  the 
Injector  (a)  Top  Left:  Axial  Velocity,  (b)  Top  Right:  Tangential  Velocity,  (c)  Bottom: 

Temperature 


Figure  104.  Comparison  of  3D  RANS  with  Experimental  Data  at  20  mm  Downstream  of  the 
Injector  (a)  Top  Left:  Axial  Velocity,  (b)  Top  Right:  Tangential  Velocity,  (c)  Bottom: 

Temperature 
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Figure  105.  Comparison  of  3D  RANS  with  Experimental  Data  at  30  mm  Downstream  of  the 
Injector  for  Velocities  and  40mm  Downstream  for  Temperature,  (a)  Top  Left:  Axial  Velocity,  (b) 
Top  Right:  Tangential  Velocity,  (c)  Bottom:  Temperature 
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Figure  106.  Comparison  of  3D  RANS  with  Experimental  Data  at  70  mm  Downstream  of  the 
Injector  for  Velocities  and  60mm  Downstream  for  Temperature,  (a)  Top  Left:  Axial  Velocity,  (b) 
Top  Right:  Tangential  Velocity,  (c)  Bottom:  Temperature 
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Figure  107.  Snapshot  of  TEC  FLAM  3D  Reacting  Flow  RANS  Simulation 


LES  is  known  to  capture  the  swirling  flow  effects  better  than  RANS.  A  LES  simulation  was 
started  in  this  phase  of  the  program  with  the  flow  modeled  just  downstream  of  the  swirl  vane 
exit.  The  flow  has  not  progressed  far  enough  to  start  averaging  and  only  instantaneous  results 
were  available.  The  better  prediction  of  the  swirl  velocity  and  the  temperature  field  is  already 
evident  from  the  results  shown  below.  LES  captures  the  much  steeper  gradients  in  the  tangential 
velocity  much  better  than  RANS.  Since  the  simulation  is  still  in  the  initial  stages,  temperatures 
near  the  walls  are  still  low  in  the  LES  result.  Hot  products  are  starting  to  recirculate  upstream, 
but  have  not  reached  the  front  of  the  combustor  on  the  centerline  and  in  the  outer  recirculation 
zone.  Plotted  results  are  shown  in  Figure  108  through  111,  while  plots  of  velocity  and 
temperature  are  shown  in  Figure  112. 
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Figure  108.  Comparison  of  3D  LES  with  Experimental  Data  at  10  mm  Downstream  of  the 
Injector,  (a)  Top  Left:  Axial  Velocity,  (b)  Top  Right:  Tangential  Velocity,  (c)  Bottom: 

Temperature 
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Figure  109.  Comparison  of  3D  LES  with  Experimental  Data  at  20  mm  Downstream  of  the 
Injector,  (a)  Top  Left:  Axial  Velocity,  (b)  Top  Right:  Tangential  Velocity,  (c)  Bottom: 

Temperature 
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Figure  110.  Comparison  of  3D  LES  with  Experimental  Data  at  30  mm  Downstream  of  the 
Injector  for  Velocities  and  40mm  Downstream  for  Temperature,  (a)  Top  Left:  Axial  Velocity,  (b) 
Top  Right:  Tangential  Velocity,  (c)  Bottom:  Temperature 


Figure  111.  Comparison  of  3D  LES  with  Experimental  Data  at  70  mm  Downstream  of  the 
Injector  for  Velocities  and  60mm  Downstream  for  Temperature,  (a)  Top  Left:  Axial  Velocity,  (b) 
Top  Right:  Tangential  Velocity,  (c)  Bottom:  Temperature 
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Figure  112.  Snapshot  of  Initial  TECFLAM  3D  Reacting  Flow  LES  Simulation 


The  extremely  large  volume  of  the  combustion  chamber  resulted  in  very  long  computational 
times.  The  results  shown  above  took  nearly  3  months  to  run.  Therefore,  it  was  decided  to 
postpone  work  on  this  case  and  focus  on  the  other  three  selected  validation  cases. 
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3.7  JSF  Combustor  Simulations 


During  the  Phase  II  SBIR  program,  RANS  (Reynolds  Averaged  Navier  Stokes)  simulations  of 
the  JSF  combustor  geometry  using  both  single-step  and  multi-step  chemistry  were  performed. 
An  LES  simulation  was  also  completed  using  single-step  kinetics.  An  LES  simulation  with 
multi-step  kinetics  was  started,  but  due  to  the  slow  convergence  behavior  of  the  solver,  was  not 
completed  during  the  Phase  II  effort.  Given  the  problem  identified  with  the  point  solver  solution 
approach  in  CFD-ACE+,  run  times  with  multi-step  kinetics  are  long.  A  single  RANS  simulation 
requires  nearly  3  weeks  of  run  time  on  16  new  computers,  while  a  single  LES  run  would  take 
nearly  4  months.  This,  and  solving  problems  related  to  the  use  of  spray  with  multi-step  kinetics 
severely  limited  the  amount  of  time  available  to  perform  JSF  combustor  simulations. 

Rolls  Royce  provided  the  combustor  geometry,  along  with  exit  profiles  of  normalized  fuel/air 
ratio.  Emissions  data  at  the  rig  conditions  simulated  for  CO  and  NOx  were  also  provided.  LES 
results  using  single-step  kinetics  resulted  in  improved  mixing  predictions  and  matched  the  NOx 
data  well.  CO  predictions  were  high  by  about  an  order  of  magnitude.  The  exit  profile 
predictions  of  fuel/air  ratio  were  also  improved  using  LES  and  matched  the  data  well.  RANS 
results  with  multi-step  kinetics  were  very  similar  to  those  with  single-step  kinetics.  However, 
CO  predictions  at  the  exit  were  essentially  equilibrium,  nearly  two  orders  of  magnitude  less  than 
the  measured  results.  The  reasons  for  this  discrepancy  are  not  clear,  but  could  be  related  to 
problems  with  the  3-step  mechanism  used.  The  PFR  kinetics  simulations  used  to  generate  the 
rate  profiles  may  not  have  been  the  best  choice.  It  is  also  unknown  how  well  the  detailed 
mechanism  used  characterized  CO  oxidation  at  the  high-pressure  conditions  of  the  combustor, 
since  validation  data  for  jet  fuels  at  high-pressure  conditions  is  extremely  limited. 

Since  the  details  of  the  JSF  combustor  simulations  are  export  controlled  and  Rolls  Royce 
proprietary,  they  have  not  been  included  in  this  version  of  the  final  report.  A  separate  export 
controlled  document  containing  this  data  was  provided  to  the  project  technical  monitor  along 
with  this  final  report. 

4.  CONCLUSIONS 

The  development  of  the  multi-step  PDF  model  represents  a  significant  improvement  in  the  state- 
of-the-art  for  practical  turbulent  combustion  modeling.  The  model  is  generally  applicable  to 
multi-step  global  mechanisms  and  results  in  a  minimal  increase  in  the  overall  computational 
expense.  Only  two  additional  transport  equations  must  be  solved  for,  regardless  of  the  number 
of  species  of  steps  in  the  mechanism.  The  model  was  successfully  tested  on  a  wide  range  of 
validation  cases.  Both  lean  premixed  and  diffusion  flame  cases  were  tested  to  showcase  the  wide 
applicability  of  the  newly  developed  turbulent  combustion  model.  Validation  cases  included 
bluff  body  stabilized  premixed,  jet  diffusion,  and  swirl  stabilized  diffusion  flames.  CFD 
simulations  were  also  performed  of  the  Rolls  Royce  JSF  combustor.  Along  with  the  new 
combustion  model,  several  other  key  improvements  were  made  to  the  LES  capabilities  in  CFD- 
ACE+.  A  dynamic  coefficient  approach  (LDVM)  was  developed  for  the  progress  and  mixture 
fraction  variance  equations.  Boundary  conditions  that  introduce  fluctuations  with  a  prescribed 
length  scale  were  implemented  to  get  more  realistic  fluctuations  on  inlet  boundaries  for  LES 
calculations. 
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Although  the  PDF  model  was  successfully  developed  and  tested,  improvements  in  the  solution  of 
the  transport  equations  will  be  essential  for  the  model  to  find  widespread  application.  During  the 
development  and  implementation  of  the  combustion  model  in  CFD-ACE+,  a  numerical  problem 
was  encountered  that  currently  limits  the  applicability  of  the  new  model.  This  is  especially  true 
for  transient  Large  Eddy  Simulation  (LES).  The  existing  method  for  solving  the  species 
transport  equations  utilizes  a  cell-by-cell  point  solver.  Specifically,  the  solver  ‘sweeps’  over  all 
of  the  computational  cells  in  the  domain  solving  the  coupled  species  equations  at  each  cell.  In  a 
sense,  this  method  holds  the  solution  values  at  all  other  cells  constant  while  solving  for  the 
species  fractions  at  a  given  cell.  While  adequate  for  small  problems,  large  problems  require  a 
large  number  of  iterations  in  order  to  assure  global  convergence.  Steady-state  problems  can 
usually  afford  a  large  number  of  global  iterations  to  achieve  ample  convergence.  However, 
transient  cases  such  as  LES  need  to  keep  the  number  of  global  iteration  per  time  step  to  a 
minimum  if  practical  run  times  are  to  be  realized. 

The  other  variables  in  the  CFD  calculations  are  solved  for  sequentially  using  a  field  solution 
procedure.  This  method  propagates  boundary  conditions  through  the  solution  domain  very 
efficiently  and  is  the  standard  approach  for  most  CFD  calculations.  Enthalpy  is  solved  for  in  this 
fashion  and  significant  problems  have  been  seen  with  the  coupling  of  the  species  (solved  for 
using  the  point  solver)  and  enthalpy.  Adopting  this  same  strategy  for  the  species  solution  is  a 
plausible  approach  to  achieve  good  coupling  of  the  enthalpy  and  species  and  was  tested.  In  this 
approach,  a  given  species  was  solved  over  the  entire  computational  domain  while  holding  all  the 
other  species  constant.  The  remaining  species  were  then  solved  for  in  turn  in  a  likewise  manner. 
This  procedure  has  the  benefit  of  fast  propagation  of  boundary  information.  However,  species 
are  highly  coupled  and  exhibit  nonlinear  behavior  and  solutions  using  this  methodology  rarely 
converge. 

The  solution  to  the  problem  would  be  to  combine  the  point  solver  and  the  field  solver  methods. 
Instead  of  solving  for  species  1  while  holding  species  2,  3,  4,  etc.  constant,  all  of  the  species 
would  be  solved  simultaneously  over  the  entire  field.  In  the  CFD-ACE+  solver,  this  would 
necessitate  that  the  field  solver  handle  more  than  current  limitation  of  solving  for  only  one 
dependent  variable  at  a  time.  Using  such  a  coupled  field  solution,  better  coupling  to  other  flow 
variables,  such  as  enthalpy,  is  possible.  This  occurs  because  the  species  information  will 
propagate  at  the  same  rate  as  the  enthalpy  information.  Also,  the  species  will  remain  implicitly 
coupled  at  each  cell.  This  solution  procedure  will  implicitly  will  require  more  CPU  time  per 
iteration,  but  will  converge  in  a  fewer  number  of  iterations.  However,  the  implementation  and 
development  of  this  capability  is  beyond  the  scope  of  the  present  program. 

In  addition  to  the  development  of  the  new  multi-step  PDF  combustion  model,  significant  work 
was  done  on  developing  the  capability  to  generate  the  needed  multi-step  mechanisms.  Reduced 
global  mechanisms  with  Arrhenius  rates,  suitable  for  CFD  calculations,  were  generated  from 
detailed  mechanisms  using  a  technique  typically  referred  to  as  Chemical  Reactor  Modeling 
(CRM).  CRM  can  be  an  effective  method  of  determining  global  Arrhenius  rate  parameters. 
Several  examples  were  shown,  demonstrating  the  flexibility  in  this  approach  with  regard  to  the 
number  of  species  and  the  number  of  steps  to  be  included  in  the  global  mechanism.  The  global 
mechanism  can  be  fine  tuned  to  combustor  operating  conditions  including  inlet  temperature, 
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equivalence  ratio  range,  and  operating  pressure.  Also,  by  using  multi-step  global  mechanisms, 
multi-step  processes  can  be  captured  that  can't  be  with  single  step  kinetics.  For  example,  multi- 
step  kinetics  can  more  accurately  capture  the  conversion  of  CH4  to  CO  and  then  to  CO2. 

There  still  exist  some  drawbacks  to  the  current  approach  for  determining  global  rate  parameters. 
For  diffusion  flame  models,  implementation  of  global  mechanisms  determined  from  CRM  can  be 
tricky  because  rate  parameters  were  not  curve  fit  for  the  entire  range  of  equivalence  ratios 
encountered.  Extrapolation  outside  of  the  bounds  where  the  rate  parameters  were  fit,  or 
application  of  the  rates  at  conditions  not  considered  in  generated  in  the  global  rates  can  result  in 
significant  errors.  Also,  the  species  profiles  should  be  included  in  the  optimization  process 
instead  of  solely  relying  on  the  temperature  profile  or  the  species  net  production  rates  (which 
have  not  proven  very  effective).  Weighting  factors  should  be  applied  to  allow  the  user  to  specify 
different  emphasis  on  temperature  or  certain  species  profiles.  Although  mechanism  generation 
can  easily  be  done  with  the  tools  developed,  generating  a  mechanism  that  behaves  realistically 
over  a  broad  range  of  conditions  and  that  accurately  mimics  the  broad  range  of  flame  properties 
necessary  (ignition  delay,  extinction  limits,  etc.)  has  proven  very  challenging  and  significant 
research  in  this  area  is  still  needed. 
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Design  of  Sample  Probes 
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Task  4. 

Experimental  Tests 
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Evaluate  Experimental  Data 
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Milestones: 

1 .  Design  of  Probes  Completed  4.  Experimental  Test  Completed 

2.  Fabrication  of  Probes  Completed  5.  Data  Evaluation  Completed 

3.  Test  Facility  Setup  Completed 


